Merge pull request #934 from YuLiu98/develop

fix(MD): init matrix in write_electric_pot function

Author: YuLiu98

source/module_md/MD_func.cpp

@@ -327,9 +327,9 @@ void MD_func::outStress(const ModuleBase::matrix &virial, const ModuleBase::matr
 }
 
 void MD_func::MDdump(const int &step, 
-		const UnitCell_pseudo &unit_in,
-		const ModuleBase::matrix &virial, 
-		const ModuleBase::Vector3<double> *force)
+        const UnitCell_pseudo &unit_in,
+        const ModuleBase::matrix &virial, 
+        const ModuleBase::Vector3<double> *force)
 {
     if(GlobalV::MY_RANK) return;
 
@@ -345,48 +345,48 @@ void MD_func::MDdump(const int &step,
         ofs.open(file.str(), ios::app);
     }
 
-	const double unit_virial = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI,3) * 1.0e-8;
-	const double unit_force = ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU;
+    const double unit_virial = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI,3) * 1.0e-8;
+    const double unit_force = ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU;
 
-	ofs << "MDSTEP:  " << step << std::endl;
-	ofs << std::setprecision(12) << std::setiosflags(ios::fixed);
+    ofs << "MDSTEP:  " << step << std::endl;
+    ofs << std::setprecision(12) << std::setiosflags(ios::fixed);
 
-	ofs << "LATTICE_CONSTANT: " << unit_in.lat0 << std::endl;
+    ofs << "LATTICE_CONSTANT: " << unit_in.lat0 << std::endl;
 
-	ofs << "LATTICE_VECTORS" << std::endl;
-	ofs << std::setw(18) << unit_in.latvec.e11 << std::setw(18) << unit_in.latvec.e12 << std::setw(18) << unit_in.latvec.e13 << std::endl; 
-	ofs << std::setw(18) << unit_in.latvec.e21 << std::setw(18) << unit_in.latvec.e22 << std::setw(18) << unit_in.latvec.e23 << std::endl;
-	ofs << std::setw(18) << unit_in.latvec.e31 << std::setw(18) << unit_in.latvec.e32 << std::setw(18) << unit_in.latvec.e33 << std::endl;
+    ofs << "LATTICE_VECTORS" << std::endl;
+    ofs << "  " << unit_in.latvec.e11 << "  " << unit_in.latvec.e12 << "  " << unit_in.latvec.e13 << std::endl; 
+    ofs << "  " << unit_in.latvec.e21 << "  " << unit_in.latvec.e22 << "  " << unit_in.latvec.e23 << std::endl;
+    ofs << "  " << unit_in.latvec.e31 << "  " << unit_in.latvec.e32 << "  " << unit_in.latvec.e33 << std::endl;
 
-	ofs << "VIRIAL (KBAR)" << std::endl;
-	for(int i=0; i<3; ++i)
-	{
-		ofs << std::setw(18) << virial(i, 0) * unit_virial 
-			<< std::setw(18) << virial(i, 1) * unit_virial 
-			<< std::setw(18) << virial(i, 2) * unit_virial << std::endl;
-	}
+    ofs << "VIRIAL (KBAR)" << std::endl;
+    for(int i=0; i<3; ++i)
+    {
+        ofs << "  " << virial(i, 0) * unit_virial 
+            << "  " << virial(i, 1) * unit_virial 
+            << "  " << virial(i, 2) * unit_virial << std::endl;
+    }
 
-	ofs << "INDEX    LABEL    POSITIONS    FORCE (eV/Angstrom)" << std::endl;
-	int index = 0;
-	for(int it=0; it<unit_in.ntype; ++it)
+    ofs << "INDEX    LABEL    POSITIONS    FORCE (eV/Angstrom)" << std::endl;
+    int index = 0;
+    for(int it=0; it<unit_in.ntype; ++it)
     {
         for(int ia=0; ia<unit_in.atoms[it].na; ++ia)
-	    {	
-		    ofs << std::setw(4) << index
-			<< std::setw(4) << unit_in.atom_label[it]
-			<< std::setw(18) << unit_in.atoms[it].tau[ia].x
-			<< std::setw(18) << unit_in.atoms[it].tau[ia].y
-			<< std::setw(18) << unit_in.atoms[it].tau[ia].z
-			<< std::setw(18) << force[index].x * unit_force 
-			<< std::setw(18) << force[index].y * unit_force 
-			<< std::setw(18) << force[index].z * unit_force << std::endl;
+        {
+            ofs << std::setw(6) << index
+            << std::setw(4) << unit_in.atom_label[it]
+            << "  " << unit_in.atoms[it].tau[ia].x
+            << "  " << unit_in.atoms[it].tau[ia].y
+            << "  " << unit_in.atoms[it].tau[ia].z
+            << "  " << force[index].x * unit_force 
+            << "  " << force[index].y * unit_force 
+            << "  " << force[index].z * unit_force << std::endl;
             index++;
-	    }
+        }
     }
 
-	ofs << std::endl;
-	ofs << std::endl;
-	ofs.close();
+    ofs << std::endl;
+    ofs << std::endl;
+    ofs.close();
 }
 
 void MD_func::getMassMbl(const UnitCell_pseudo &unit_in, 

source/src_io/write_pot.cpp

@@ -268,7 +268,12 @@ void Potential::write_elecstat_pot(const std::string &fn, const std::string &fn_
     //==========================================
     for (int ir = 0;ir < GlobalC::pw.nrxx;ir++)
     {
-        v_elecstat[ir] = Porter[ir].real() + this->vltot[ir] + v_efield(0, ir);
+        v_elecstat[ir] = Porter[ir].real() + this->vltot[ir];
+
+        if (GlobalV::EFIELD && GlobalV::DIPOLE)
+        {
+            v_elecstat[ir] += v_efield(0, ir);
+        }
     }
 
     //-------------------------------------------

source/src_pw/forces.cpp

@@ -76,6 +76,10 @@ void Forces::init(ModuleBase::matrix& force)
     {
         force_e.create(GlobalC::ucell.nat, 3);
         Efield::compute_force(GlobalC::ucell, force_e);
+        if(GlobalV::TEST_FORCE)
+        {
+            Forces::print("EFIELD      FORCE (Ry/Bohr)", force_e);
+        }
     }
 
     //impose total force = 0