Merge update_MD branch to develop branch

doc/input-main.md

@@ -43,7 +43,7 @@
 
     - [Molecular dynamics](#molecular-dynamics)
 
-        [md_type](#md-type) | [md_potential](#md-potential) | [md_rstmd](#md-rstmd) | [md_dt](#md_dt) | [md_t](#md-t) | [md_qmass](#md-qmass) | [md_dumpmdfred](#md-dumpmdfred) | [md_fixtemperature](#md-fixtemperature) | [NVT_control](#nvt-control) | [NVT_tau](#nvt-tau) | [MNHC](#mnhc) | [md_ediff](#md-ediff) | [md_ediffg](#md-ediffg) | [rcut_lj](#rcut_lj) | [epsilon_lj](#epsilon_lj) | [sigma_lj](#sigma_lj)
+        [md_type](#md-type) | [md_potential](#md-potential) | [md_rstmd](#md-rstmd) | [md_dt](#md_dt) | [md_t](#md-t) | [md_qmass](#md-qmass) | [md_dumpmdfred](#md-dumpmdfred) | [md_tfreq](#md-tfreq) | [md_fixtemperature](#md-fixtemperature) | [NVT_control](#nvt-control) | [NVT_tau](#nvt-tau) | [MNHC](#mnhc) | [md_ediff](#md-ediff) | [md_ediffg](#md-ediffg) | [rcut_lj](#rcut_lj) | [epsilon_lj](#epsilon_lj) | [sigma_lj](#sigma_lj) | [direction](#direction) | [velocity](#velocity) | [viscosity](#viscosity) | [tscale](#tscale) | [damp](#damp)
 
     - [DFT+U correction](#DFT_U-correction)
 
@@ -1111,10 +1111,13 @@ This part of variables are used to control the molecular dynamics calculations.
 - md_type<a id="md-type"></a>
     - *Type*: Integer
     - *Description*: control the ensemble to run md.
-        - 0: When set to 0, ABACUS will use NVE ensemble;
-        - 1: When set to 1, ABACUS will use NVT ensemble with Nose Hoover method;
-        - 2: When set to 2, ABACUS will use NVT ensemble with Velosity Scaling method;
-    - *Default*: 1
+        - -1: FIRE method to relax;
+        - 0: NVE ensemble;
+        - 1: NVT ensemble with Anderson thermostat;
+        - 2: NVT ensemble with Nose Hoover Chain;
+        - 3: NVT ensemble with Langevin method;
+        - 4: MSST method; 
+    - *Default*: 2
 
     [back to top](#input-file)
 
@@ -1131,15 +1134,15 @@ This part of variables are used to control the molecular dynamics calculations.
 - md_rstmd<a id="md-rstmd"></a>
     - *Type*: Boolean
     - *Description*: to control whether restart md.
-        - 0:When set to 0, ABACUS will calculate md normolly.
-        - 1:When set to 1, ABACUS will calculate md from last step in your test before.
+        - 0: When set to 0, ABACUS will calculate md normolly.
+        - 1: When set to 1, ABACUS will calculate md from last step in your test before.
     - *Default*: 0
 
     [back to top](#input-file)
 - md_dt<a id="md_dt"></a>
     - *Type*: Double
     - *Description*: This is the time step(fs) used in md simulation .
-    - *Default*: No default
+    - *Default*: 1
 
     [back to top](#input-file)
 - md_tfirst & md_tlast<a id="md-t"></a>
@@ -1150,8 +1153,8 @@ This part of variables are used to control the molecular dynamics calculations.
     [back to top](#input-file)
 - md_qmass<a id="md-qmass"></a>
     - *Type*: Double
-    - *Description*: Inertia of extended system variable. Used only when md_type is 1 or 2, you should set a number which is larger than 0. If you want to autoset this by ABACUS,just set it to 0.
-    - *Default*: 0
+    - *Description*: Inertia of extended system variable. Used only when md_type is 4, you should set a number which is larger than 0. Note that Qmass of NHC is set by md_tfreq.
+    - *Default*: No default
 
     [back to top](#input-file)
 - md_dumpmdfred<a id="md-dumpmdfred"></a>
@@ -1161,6 +1164,14 @@ This part of variables are used to control the molecular dynamics calculations.
 
     [back to top](#input-file)
 
+- md_tfreq<a id="md-tfreq"></a>
+    - *Type*: Real
+    - *Description*: 
+        - Oscillation frequency, used to determine Qmass of NHC; 
+        - 1/(md_tfreq*md_dt) is collision probability in Anderson method.
+    - *Default*: 1.0
+
+    [back to top](#input-file)
 
 - md_fixtemperature<a id="md-fixtemperature"></a>
     - *Type*: Integer
@@ -1230,6 +1241,41 @@ This part of variables are used to control the molecular dynamics calculations.
 
     [back to top](#input-file)
 
+- direction<a id="direction"></a>
+    - *Type*: Integer
+    - *Description*: the direction of shock wave for MSST.
+    - *Default*: 2 (z direction)
+
+    [back to top](#input-file)
+
+- velocity<a id="velocity"></a>
+    - *Type*: Real
+    - *Description*: the velocity of shock wave (\AA/fs) for MSST.
+    - *Default*: 0
+
+    [back to top](#input-file)
+
+- viscosity<a id="viscosity"></a>
+    - *Type*: Real
+    - *Description*: artificial viscosity (mass/length/time) for MSST.
+    - *Default*: 0
+
+    [back to top](#input-file)
+
+- tscale<a id="tscale"></a>
+    - *Type*: Real
+    - *Description*: reduction in initial temperature (0~1) used to compress volume in MSST.
+    - *Default*: 0
+
+    [back to top](#input-file)
+
+- damp<a id="damp"></a>
+    - *Type*: Real
+    - *Description*: damping parameter (fs) used to add force in Langevin method.
+    - *Default*: 1
+
+    [back to top](#input-file)
+
 ### DFT+U correction
 This part of variables are used to control DFT+U correlated parameters
 - dft_plus_u<a id="dft_plus_u"></a>

source/Makefile.Objects

@@ -55,9 +55,6 @@ chi0_hilbert.o\
 chi0_standard.o\
 epsilon0_pwscf.o\
 epsilon0_vasp.o\
-MD_basic.o\
-MD_thermo.o\
-MD_fire.o\
 MD_func.o\
 exx_lip.o\
 soc.o\
@@ -298,6 +295,13 @@ run_md_classic.o\
 LJ_potential.o\
 DP_potential.o\
 cmd_neighbor.o\
+verlet.o\
+NVE.o\
+MSST.o\
+NVT_ADS.o\
+NVT_NHC.o\
+FIRE.o\
+Langevin.o\
 dos.o \
 inverse_matrix.o \
 energy.o \

source/input.cpp

@@ -1218,6 +1218,30 @@ bool Input::Read(const std::string &fn)
 		{
 			read_value(ifs, mdp.md_potential);
 		}
+		else if (strcmp("direction",word) == 0)
+		{
+			read_value(ifs, mdp.direction);
+		}
+		else if (strcmp("velocity",word) == 0)
+		{
+			read_value(ifs, mdp.velocity);
+		}
+		else if (strcmp("viscosity",word) == 0)
+		{
+			read_value(ifs, mdp.viscosity);
+		}
+		else if (strcmp("tscale",word) == 0)
+		{
+			read_value(ifs, mdp.tscale);
+		}
+		else if (strcmp("md_tfreq",word) == 0)
+		{
+			read_value(ifs, mdp.tfreq);
+		}
+		else if (strcmp("md_damp",word) == 0)
+		{
+			read_value(ifs, mdp.damp);
+		}
 //----------------------------------------------------------
 // tddft
 // Fuxiang He add 2016-10-26
@@ -2146,6 +2170,12 @@ void Input::Bcast()
 		Parallel_Common::bcast_double(mdp.epsilon_lj);
 		Parallel_Common::bcast_double(mdp.sigma_lj);
 		Parallel_Common::bcast_string(mdp.md_potential);
+		Parallel_Common::bcast_int(mdp.direction);
+		Parallel_Common::bcast_double(mdp.velocity);
+		Parallel_Common::bcast_double(mdp.viscosity);
+		Parallel_Common::bcast_double(mdp.tscale);
+		Parallel_Common::bcast_double(mdp.tfreq);
+		Parallel_Common::bcast_double(mdp.damp);
 /* 	// Peize Lin add 2014-04-07
 	Parallel_Common::bcast_bool( vdwD2 );
 	Parallel_Common::bcast_double( vdwD2_scaling );

source/module_base/constants.h

@@ -72,7 +72,7 @@ const double ANGSTROM_AU      = 1.8897270;		// au
 //const double AU_TERAHERTZ     = 2.418e-5;		// THz
 //const double TERAHERTZ        = 2.418e-5;		// from au to THz
 //const double AU_SEC           = 2.4189e-17;		// sec
-//const double AU_PS            = 2.4189e-5;		// sec
+const double AU_to_FS           = 2.418884326505e-2;		// from a.u. to fs
 //const double rhothr = 1.0e-5;					// tolerance
 //const double gsmall = 1.0e-12;
 const double e2 = 2.0;							// the square of the electron charge

source/module_base/matrix3.h

@@ -15,7 +15,7 @@ class Matrix3
 {
 	/* data */
 public:
-	// element eij:i_column,j_row
+	// element eij: i_row, j_column
 	double e11, e12, e13, e21, e22, e23, e31, e32, e33;
 
 	/* Constructors and destructor */

source/module_base/vector3.h

@@ -32,6 +32,9 @@ class Vector3
 	Vector3<T>& operator/=(const T &s)          { x/=s; y/=s; z/=s;        return *this; }
 	Vector3<T>  operator -() const              { return Vector3<T>(-x,-y,-z); }			// Peize Lin add 2017-01-10
 
+	T operator[](int index)const { return (&x)[index]; }
+	T& operator[](int index) { return (&x)[index]; }
+
 	T norm2(void) const	{ return x*x + y*y + z*z; }
 	T norm(void) const	{ return sqrt(norm2()); }
 	Vector3<T>& normalize(void){ const T m=norm(); x/=m; y/=m; z/=m; return *this; }	// Peize Lin update return 2019-09-08

source/module_cell/atom_spec.cpp

@@ -78,7 +78,7 @@ void Atom::set_index(void)
     return;
 }
 
-void Atom::print_Atom(std::ofstream &ofs, output &outp)
+void Atom::print_Atom(std::ofstream &ofs)
 {
     //OUT(ofs,"print_Atom()");
     ModuleBase::GlobalFunc::OUT(ofs,"label",label);
@@ -95,7 +95,7 @@ void Atom::print_Atom(std::ofstream &ofs, output &outp)
     //===================
     this->print_atom(ofs);
 
-    outp.printv31_d(ofs,"atom_position(cartesian)",tau,na);
+    output::printv31_d(ofs,"atom_position(cartesian)",tau,na);
     /*
     for (int i = 0;i < na;i++)
     {
@@ -134,11 +134,11 @@ void Atom::bcast_atom(void)
         assert(na!=0);
         delete[] tau;
 		delete[] taud;
-	delete[] vel;
+	    delete[] vel;
         delete[] mag;
         tau = new ModuleBase::Vector3<double>[na];
 		taud = new ModuleBase::Vector3<double>[na];
-	vel = new ModuleBase::Vector3<double>[na];
+	    vel = new ModuleBase::Vector3<double>[na];
         mag = new double[na];
         angle1 = new double[na];
         angle2 = new double[na];
@@ -153,9 +153,9 @@ void Atom::bcast_atom(void)
         Parallel_Common::bcast_double( taud[i].x );
         Parallel_Common::bcast_double( taud[i].y );
         Parallel_Common::bcast_double( taud[i].z );
-	Parallel_Common::bcast_double( vel[i].x );
-	Parallel_Common::bcast_double( vel[i].y );
-	Parallel_Common::bcast_double( vel[i].z );
+	    Parallel_Common::bcast_double( vel[i].x );
+	    Parallel_Common::bcast_double( vel[i].y );
+	    Parallel_Common::bcast_double( vel[i].z );
         Parallel_Common::bcast_double( mag[i] );
         Parallel_Common::bcast_double(angle1[i]);
         Parallel_Common::bcast_double(angle2[i]);

source/module_cell/atom_spec.h

@@ -41,7 +41,7 @@ class Atom: public Atom_pseudo
     ModuleBase::Vector3<double> *m_loc_;
 
 
-    void print_Atom(std::ofstream &ofs, output &outp);
+    void print_Atom(std::ofstream &ofs);
 #ifdef __MPI
     void bcast_atom(void);
     void bcast_atom2(void);

source/module_cell/pseudo_nc.cpp

@@ -99,15 +99,15 @@ void pseudo_nc::set_pseudo_nc(const Pseudopot_upf &upf)
 } // end subroutine set_pseudo_upf
 
 
-void pseudo_nc::print_pseudo_nc(std::ofstream &ofs, output &outp)
+void pseudo_nc::print_pseudo_nc(std::ofstream &ofs)
 {
-	print_pseudo_vl(ofs, outp);
+	print_pseudo_vl(ofs);
 	ofs << "\n pseudo_nc : ";
 	ofs << "\n kkbeta	" << kkbeta;
 	ofs << "\n nh  " << nh;
-	outp.printr1_d(ofs, " lll : ", lll, nbeta);
-	outp.printrm(ofs, " betar : ", betar);
-	outp.printrm(ofs, " dion : ", dion);
+	output::printr1_d(ofs, " lll : ", lll, nbeta);
+	output::printrm(ofs, " betar : ", betar);
+	output::printrm(ofs, " dion : ", dion);
 	ofs << "\n ----------------------";
 }
 
@@ -345,31 +345,31 @@ void pseudo_nc::set_pseudo_vl(const Pseudopot_upf &upf)
 } 
 
 
-void pseudo_nc::print_pseudo_atom(std::ofstream &ofs, output &outp)
+void pseudo_nc::print_pseudo_atom(std::ofstream &ofs)
 {
-	print_pseudo_h(ofs, outp);
+	print_pseudo_h(ofs);
 	ofs << "\n pseudo_atom : ";
 	ofs << "\n msh	" << msh;
 //	ofs	<< "\n nchi	" << nchi;
-	outp.printr1_d(ofs, " r : ", r, mesh);
-	outp.printr1_d(ofs, " rab : ", rab, mesh);
-	outp.printr1_d(ofs, " rho_atc : ", rho_atc, mesh);
-	outp.printr1_d(ofs, " rho_at : ", rho_at, mesh);
-	outp.printr1_d(ofs," jchi : ", jchi, nchi);
-	outp.printrm(ofs, " chi : ", chi);
+	output::printr1_d(ofs, " r : ", r, mesh);
+	output::printr1_d(ofs, " rab : ", rab, mesh);
+	output::printr1_d(ofs, " rho_atc : ", rho_atc, mesh);
+	output::printr1_d(ofs, " rho_at : ", rho_at, mesh);
+	output::printr1_d(ofs," jchi : ", jchi, nchi);
+	output::printrm(ofs, " chi : ", chi);
 	ofs << "\n ----------------------";
 }
 
 
-void pseudo_nc::print_pseudo_vl(std::ofstream &ofs, output &outp)
+void pseudo_nc::print_pseudo_vl(std::ofstream &ofs)
 {
 	ofs << "\n pseudo_vl:";
-	print_pseudo_atom(ofs, outp);
-	outp.printr1_d(ofs, "vloc_at : ", vloc_at, mesh);
+	print_pseudo_atom(ofs);
+	output::printr1_d(ofs, "vloc_at : ", vloc_at, mesh);
 	ofs << "\n ----------------------------------- ";
 }
 
-void pseudo_nc::print_pseudo_h(std::ofstream &ofs, output &outp)
+void pseudo_nc::print_pseudo_h(std::ofstream &ofs)
 {
     ofs << "\n pseudo_info :";
     ofs << "\n nv       " << nv;
@@ -387,8 +387,8 @@ void pseudo_nc::print_pseudo_h(std::ofstream &ofs, output &outp)
     ofs << "\n nchi " << nchi;
     ofs << "\n nbeta    " << nbeta;
 //  out.printr1_d(ofs," els: ", els, nchi);
-    outp.printr1_d(ofs, " lchi: ", lchi, nchi);
-    outp.printr1_d(ofs, " oc: ", oc, nchi);
+    output::printr1_d(ofs, " lchi: ", lchi, nchi);
+    output::printr1_d(ofs, " oc: ", oc, nchi);
     ofs << "\n ----------------------";
 }
 

source/module_cell/pseudo_nc.h

@@ -76,10 +76,10 @@ class pseudo_nc
 	void set_pseudo_vl(const Pseudopot_upf &upf);
 	void set_pseudo_nc(const Pseudopot_upf &upf);
 
-	void print_pseudo_h(std::ofstream &ofs, output &outp);
-	void print_pseudo_atom(std::ofstream &ofs, output &outp);
-	void print_pseudo_vl(std::ofstream &ofs, output &outp);
-	void print_pseudo_nc(std::ofstream &ofs, output &outp);
+	void print_pseudo_h(std::ofstream &ofs);
+	void print_pseudo_atom(std::ofstream &ofs);
+	void print_pseudo_vl(std::ofstream &ofs);
+	void print_pseudo_nc(std::ofstream &ofs);
 
 };