Skip to content

Commit 4b9442b

Browse files
YuLiu98YuLiu98
authored
Merge pull request #966 from YuLiu98/develop
fix : fix md bug when cal_stress=0
2 parents 01e8920 + 7dd4910 commit 4b9442b

File tree

20 files changed

+99
-66
lines changed

20 files changed

+99
-66
lines changed

‎source/module_md/FIRE.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -86,9 +86,9 @@ void FIRE::second_half()
8686
ModuleBase::timer::tick("FIRE", "second_half");
8787
}
8888

89-
void FIRE::outputMD(std::ofstream &ofs)
89+
void FIRE::outputMD(std::ofstream &ofs, bool cal_stress)
9090
{
91-
Verlet::outputMD(ofs);
91+
Verlet::outputMD(ofs, cal_stress);
9292

9393
ofs << " LARGEST GRAD (eV/A) : "
9494
<< max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl;

‎source/module_md/FIRE.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ class FIRE : public Verlet
1212
void setup(ModuleESolver::ESolver *p_ensolve);
1313
void first_half();
1414
void second_half();
15-
void outputMD(std::ofstream &ofs);
15+
void outputMD(std::ofstream &ofs, bool cal_stress);
1616
void write_restart();
1717
void restart();
1818
void check_force();

‎source/module_md/Langevin.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -45,9 +45,9 @@ void Langevin::second_half()
4545
ModuleBase::timer::tick("Langevin", "second_half");
4646
}
4747

48-
void Langevin::outputMD(std::ofstream &ofs)
48+
void Langevin::outputMD(std::ofstream &ofs, bool cal_stress)
4949
{
50-
Verlet::outputMD(ofs);
50+
Verlet::outputMD(ofs, cal_stress);
5151
}
5252

5353
void Langevin::write_restart()

‎source/module_md/Langevin.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ class Langevin : public Verlet
1212
void setup(ModuleESolver::ESolver *p_ensolve);
1313
void first_half();
1414
void second_half();
15-
void outputMD(std::ofstream &ofs);
15+
void outputMD(std::ofstream &ofs, bool cal_stress);
1616
void write_restart();
1717
void restart();
1818
void post_force();

‎source/module_md/MD_func.cpp‎

Lines changed: 40 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -75,36 +75,39 @@ void MD_func::kinetic_stress(
7575
// This function calculates the classical kinetic energy of atoms
7676
// and its contribution to stress.
7777
//----------------------------------------------------------------------------
78-
kinetic = MD_func::GetAtomKE(unit_in.nat, vel, allmass);
78+
kinetic = MD_func::GetAtomKE(unit_in.nat, vel, allmass);
7979

80-
ModuleBase::matrix temp;
81-
temp.create(3,3); // initialize
80+
if(GlobalV::CAL_STRESS)
81+
{
82+
ModuleBase::matrix temp;
83+
temp.create(3,3); // initialize
8284

83-
for(int ion=0; ion<unit_in.nat; ++ion)
84-
{
85-
for(int i=0; i<3; ++i)
86-
{
87-
for(int j=i; j<3; ++j)
88-
{
89-
temp(i, j) += allmass[ion] * vel[ion][i] * vel[ion][j];
90-
}
91-
}
92-
}
85+
for(int ion=0; ion<unit_in.nat; ++ion)
86+
{
87+
for(int i=0; i<3; ++i)
88+
{
89+
for(int j=i; j<3; ++j)
90+
{
91+
temp(i, j) += allmass[ion] * vel[ion][i] * vel[ion][j];
92+
}
93+
}
94+
}
9395

94-
for(int i=0; i<3; ++i)
95-
{
96-
for(int j=0; j<3; ++j)
97-
{
98-
if(j<i)
99-
{
100-
stress(i, j) = stress(j, i);
101-
}
102-
else
103-
{
104-
stress(i, j) = temp(i, j)/unit_in.omega;
105-
}
106-
}
107-
}
96+
for(int i=0; i<3; ++i)
97+
{
98+
for(int j=0; j<3; ++j)
99+
{
100+
if(j<i)
101+
{
102+
stress(i, j) = stress(j, i);
103+
}
104+
else
105+
{
106+
stress(i, j) = temp(i, j)/unit_in.omega;
107+
}
108+
}
109+
}
110+
}
108111
}
109112

110113
// Read Velocity from STRU liuyu 2021-09-24
@@ -239,6 +242,7 @@ void MD_func::force_virial(
239242

240243
if(mdp.md_ensolver == "LJ")
241244
{
245+
GlobalV::CAL_STRESS = 1;
242246
bool which_method = unit_in.judge_big_cell();
243247
if(which_method)
244248
{
@@ -272,6 +276,7 @@ void MD_func::force_virial(
272276
}
273277
else if(mdp.md_ensolver == "DP")
274278
{
279+
GlobalV::CAL_STRESS = 1;
275280
DP_potential::DP_pot(unit_in, potential, force, stress);
276281
}
277282
#ifndef __CMD
@@ -358,12 +363,15 @@ void MD_func::MDdump(const int &step,
358363
ofs << " " << unit_in.latvec.e21 << " " << unit_in.latvec.e22 << " " << unit_in.latvec.e23 << std::endl;
359364
ofs << " " << unit_in.latvec.e31 << " " << unit_in.latvec.e32 << " " << unit_in.latvec.e33 << std::endl;
360365

361-
ofs << "VIRIAL (KBAR)" << std::endl;
362-
for(int i=0; i<3; ++i)
366+
if(GlobalV::CAL_STRESS)
363367
{
364-
ofs << " " << virial(i, 0) * unit_virial
365-
<< " " << virial(i, 1) * unit_virial
366-
<< " " << virial(i, 2) * unit_virial << std::endl;
368+
ofs << "VIRIAL (KBAR)" << std::endl;
369+
for(int i=0; i<3; ++i)
370+
{
371+
ofs << " " << virial(i, 0) * unit_virial
372+
<< " " << virial(i, 1) * unit_virial
373+
<< " " << virial(i, 2) * unit_virial << std::endl;
374+
}
367375
}
368376

369377
ofs << "INDEX LABEL POSITIONS FORCE (eV/Angstrom)" << std::endl;

‎source/module_md/MSST.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ MSST::MSST(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para
1010
{
1111
std::cout << "MSST" << std::endl;
1212

13+
GlobalV::CAL_STRESS = 1;
1314
mdp.msst_qmass = mdp.msst_qmass / pow(ModuleBase::ANGSTROM_AU, 4) / pow(ModuleBase::AU_to_MASS, 2);
1415
mdp.msst_vel = mdp.msst_vel * ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS;
1516
mdp.msst_vis = mdp.msst_vis / ModuleBase::AU_to_MASS / ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS;
@@ -158,9 +159,9 @@ void MSST::second_half()
158159
ModuleBase::timer::tick("MSST", "second_half");
159160
}
160161

161-
void MSST::outputMD(std::ofstream &ofs)
162+
void MSST::outputMD(std::ofstream &ofs, bool cal_stress)
162163
{
163-
Verlet::outputMD(ofs);
164+
Verlet::outputMD(ofs, cal_stress);
164165
}
165166

166167
void MSST::write_restart()

‎source/module_md/MSST.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ class MSST : public Verlet
1212
void setup(ModuleESolver::ESolver *p_ensolve);
1313
void first_half();
1414
void second_half();
15-
void outputMD(std::ofstream &ofs);
15+
void outputMD(std::ofstream &ofs, bool cal_stress);
1616
void write_restart();
1717
void restart();
1818
double extra_term();

‎source/module_md/NVE.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -36,9 +36,9 @@ void NVE::second_half()
3636
ModuleBase::timer::tick("NVE", "second_half");
3737
}
3838

39-
void NVE::outputMD(std::ofstream &ofs)
39+
void NVE::outputMD(std::ofstream &ofs, bool cal_stress)
4040
{
41-
Verlet::outputMD(ofs);
41+
Verlet::outputMD(ofs, cal_stress);
4242
}
4343

4444
void NVE::write_restart()

‎source/module_md/NVE.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ class NVE : public Verlet
1212
void setup(ModuleESolver::ESolver *p_ensolve);
1313
void first_half();
1414
void second_half();
15-
void outputMD(std::ofstream &ofs);
15+
void outputMD(std::ofstream &ofs, bool cal_stress);
1616
void write_restart();
1717
void restart();
1818

‎source/module_md/NVT_ADS.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -65,9 +65,9 @@ void NVT_ADS::second_half()
6565
ModuleBase::timer::tick("NVT_ADS", "second_half");
6666
}
6767

68-
void NVT_ADS::outputMD(std::ofstream &ofs)
68+
void NVT_ADS::outputMD(std::ofstream &ofs, bool cal_stress)
6969
{
70-
Verlet::outputMD(ofs);
70+
Verlet::outputMD(ofs, cal_stress);
7171
}
7272

7373
void NVT_ADS::write_restart()

0 commit comments

Comments
 (0)