diff --git a/doc/examples/md.md b/doc/examples/md.md index 69611611e6..469e299cb9 100644 --- a/doc/examples/md.md +++ b/doc/examples/md.md @@ -34,7 +34,6 @@ mixing_beta 0.4 charge_extrap second-order md_mdtype 1 //choose ensemble -md_potential 0 //choose potential 0:Pseudopotential 1:LJ potential md_dt 1 //time step md_tfirst 700 //the first target temperature md_rstmd 0 //whether restart md @@ -47,7 +46,6 @@ These MD parameters means that ABACUS will use NVT ensemble with Nosé-hoover th Note: *Please turn off symmetry when do MD simulation.* - md_mdtype : 0, NVE; 1, NVT; 2, velocity scaling -- md_potential : 0, Pseudopotential 1, LJ potential - md_dt : time step in md simulation (fs) - md_tfirst : target temperature in md simulation(K), you should set parameter md_tlast and md_fixtemperature when you want to change temperature during md simulation. - md_rstmd : 0, no need of restart ; 1, restart with restart file, you must repalce STRU file with STRU_MD before you run the restart task. diff --git a/source/driver.cpp b/source/driver.cpp index ae02f8dd48..76a679d851 100644 --- a/source/driver.cpp +++ b/source/driver.cpp @@ -81,17 +81,8 @@ void Driver::atomic_world(void) // lcao_in_pw: LCAO expaned by plane wave basis set // lcao: linear combination of atomic orbitals //-------------------------------------------------- -#ifndef __NOMD //qianrui do not want to add md to pw module temporarily before md is compeletly built. - if(CALCULATION=="md" || CALCULATION=="md-sto") // Yu Liu 2021-07-12 - { - Run_md::md_line(); - } - else if(BASIS_TYPE=="pw" || BASIS_TYPE=="lcao_in_pw") -#else if(BASIS_TYPE=="pw" || BASIS_TYPE=="lcao_in_pw") -#endif { - Run_pw::plane_wave_line(); } #ifdef __LCAO diff --git a/source/input.cpp b/source/input.cpp index b69c171816..96321be35f 100644 --- a/source/input.cpp +++ b/source/input.cpp @@ -165,6 +165,7 @@ void Input::Default(void) search_radius=-1.0; // unit: a.u. -1.0 has no meaning. search_pbc=true; symmetry=false; + set_vel=false; mlwf_flag=false; force=0; force_set=false; @@ -710,6 +711,10 @@ bool Input::Read(const string &fn) else if (strcmp("symmetry", word) == 0) { read_value(ifs, symmetry); + } + else if (strcmp("set_vel", word) == 0) + { + read_value(ifs, set_vel); } else if (strcmp("mlwf_flag", word) == 0) { @@ -1148,10 +1153,10 @@ bool Input::Read(const string &fn) { read_value(ifs, mdp.mdtype); } - else if (strcmp("md_potential",word) == 0) - { - read_value(ifs, mdp.md_potential); - } + //else if (strcmp("md_potential",word) == 0) + //{ + // read_value(ifs, mdp.md_potential); + //} else if (strcmp("NVT_tau",word) == 0) { read_value(ifs, mdp.NVT_tau); @@ -2010,6 +2015,7 @@ void Input::Bcast() Parallel_Common::bcast_bool( search_pbc ); Parallel_Common::bcast_double( search_radius ); Parallel_Common::bcast_bool( symmetry ); + Parallel_Common::bcast_bool( set_vel ); //liuyu 2021-07-14 Parallel_Common::bcast_bool( mlwf_flag ); Parallel_Common::bcast_int( force ); Parallel_Common::bcast_bool( force_set ); @@ -2127,7 +2133,7 @@ void Input::Bcast() */ //zheng daye add 2014/5/5 Parallel_Common::bcast_int(mdp.mdtype); - Parallel_Common::bcast_int(mdp.md_potential); + //Parallel_Common::bcast_int(mdp.md_potential); Parallel_Common::bcast_double(mdp.NVT_tau); Parallel_Common::bcast_int(mdp.NVT_control); Parallel_Common::bcast_double(mdp.dt); diff --git a/source/input.h b/source/input.h index 671173f02a..8afc2616ac 100644 --- a/source/input.h +++ b/source/input.h @@ -51,6 +51,8 @@ class Input double emin_sto; string stotype; + bool set_vel; // read velocity from STRU or not liuyu 2021-07-14 + bool symmetry; // turn on symmetry or not int npool; // ecch pool is for one k point diff --git a/source/module_cell/atom_spec.cpp b/source/module_cell/atom_spec.cpp index b3210dc03e..b50cf3ac33 100644 --- a/source/module_cell/atom_spec.cpp +++ b/source/module_cell/atom_spec.cpp @@ -13,6 +13,7 @@ Atom::Atom() stapos_wf = 0; tau = new Vector3[1]; taud = new Vector3[1]; + vel = new Vector3[1]; mag = new double[1]; l_nchi = new int[1]; iw2l = new int[1]; @@ -26,6 +27,7 @@ Atom::~Atom() { delete[] tau; delete[] taud; + delete[] vel; delete[] mag; delete[] l_nchi; delete[] iw2l; diff --git a/source/module_cell/atom_spec.h b/source/module_cell/atom_spec.h index 91635fc676..1b7ae470c1 100644 --- a/source/module_cell/atom_spec.h +++ b/source/module_cell/atom_spec.h @@ -33,6 +33,7 @@ class Atom: public Atom_pseudo string label; // atomic symbol Vector3 *tau;// Cartesian coordinates of each atom in this type. Vector3 *taud;// Direct coordinates of each atom in this type. + Vector3 *vel;// velocities of each atom in this type. double* mag; double angle1;//spin angle, added by zhengdy-soc diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp index 49b15c1a13..22ddc9b77d 100644 --- a/source/module_cell/read_atoms.cpp +++ b/source/module_cell/read_atoms.cpp @@ -398,10 +398,12 @@ bool UnitCell_pseudo::read_atom_positions(ifstream &ifpos) { delete[] atoms[it].tau; delete[] atoms[it].taud; + delete[] atoms[it].vel; delete[] atoms[it].mbl; delete[] atoms[it].mag; atoms[it].tau = new Vector3[na]; atoms[it].taud = new Vector3[na]; + atoms[it].vel = new Vector3[na]; atoms[it].mbl = new Vector3[na]; atoms[it].mag = new double[na]; atoms[it].mass = this->atom_mass[it]; //mohan add 2011-11-07 @@ -436,6 +438,14 @@ bool UnitCell_pseudo::read_atom_positions(ifstream &ifpos) { ifpos >> v.x >> v.y >> v.z >> mv.x >> mv.y >> mv.z; + if(INPUT.set_vel) + { + ifpos >> atoms[it].vel[ia].x >> atoms[it].vel[ia].y >> atoms[it].vel[ia].z; + } + else + { + atoms[it].vel[ia].set(0,0,0); + } string mags; std::getline( ifpos, mags ); // change string to double. @@ -675,7 +685,8 @@ void UnitCell_pseudo::print_stru_file(const string &fn, const int &type)const for(int ia=0; ia[ucell.nat]; + vel=new Vector3[ucell.nat]; //initialize velocity of atoms from STRU Yu Liu 2021-07-14 + int iat=0; + for(int it=0; it[ucell.nat]; tauDirectChange=new Vector3[ucell.nat]; allmass=new double[ucell.nat]; diff --git a/source/module_md/MD_parameters.h b/source/module_md/MD_parameters.h index e3b5b4fc54..f0ebf5b4bd 100644 --- a/source/module_md/MD_parameters.h +++ b/source/module_md/MD_parameters.h @@ -7,7 +7,7 @@ class MD_parameters MD_parameters() { mdtype=1; - md_potential=0; + //md_potential=0; NVT_tau=0; NVT_control=1; dt=1; @@ -27,7 +27,7 @@ class MD_parameters int rstMD; // 1 : restart MD, vel.restart and ion.restart files will be read // 0 : not restart from ion and vel files int mdtype; // 1: NVT:Nose Hoover, 2:NVT:velocity scaling 3: NPT, 0: NVE - int md_potential; // 0: Pseudopotential, 1: LJ potential Yu Liu 2021-07-12 + //int md_potential; // 0: Pseudopotential, 1: LJ potential Yu Liu 2021-07-12 double Qmass; // Inertia of extended system variable double dt ; // Time increment (hbar/E_hartree) double tfirst; //Temperature (in Hartree, 1 Hartree ~ 3E5 K) diff --git a/source/module_md/run_md.cpp b/source/module_md/run_md.cpp index 3b6d84a670..006025fa1a 100644 --- a/source/module_md/run_md.cpp +++ b/source/module_md/run_md.cpp @@ -9,6 +9,7 @@ #include "../src_io/cal_test.h" #include "../src_io/print_info.h" #include "../src_pw/symmetry.h" +#include "../module_neighbor/sltk_atom_arrange.h" Run_md::Run_md(){} Run_md::~Run_md(){} @@ -19,140 +20,32 @@ void Run_md::md_line(void) TITLE("Run_md","md_line"); timer::tick("Run_md","md_line"); - if(INPUT.mdp.md_potential==0) - { - Run_md::ai_md_line(); - } - else - { Run_md::classic_md_line(); - } timer::tick("Run_md","md_line"); return; } - -void Run_md::ai_md_line(void) -{ - // Setup the unitcell. - // improvement: a) separating the first reading of the atom_card and subsequent - // cell relaxation. b) put NLOCAL and NBANDS as input parameters - ucell.setup_cell( global_pseudo_dir, out, global_atom_card, ofs_running); - - // setup NBANDS - // Yu Liu add 2021-07-03 - CHR.cal_nelec(); - - // mohan add 2010-09-06 - // Yu Liu move here 2021-06-27 - // because the number of element type - // will easily be ignored, so here - // I warn the user again for each type. - for(int it=0; itmd_allocate_ions(); + + MD_basic mdb(INPUT.mdp, ucell); + int mdtype = INPUT.mdp.mdtype; + + this->istep = 1; + bool stop = false; + + while (istep <= NSTEP && !stop) + { + time_t estart = time(NULL); + + if (OUT_LEVEL == "ie") + { + cout << " -------------------------------------------" << endl; + cout << " STEP OF MOLECULAR DYNAMICS : " << istep << endl; + cout << " -------------------------------------------" << endl; + ofs_running << " -------------------------------------------" << endl; + ofs_running << " STEP OF MOLECULAR DYNAMICS : " << istep << endl; + ofs_running << " -------------------------------------------" << endl; + } + + this->update_pos_classic(); + + time_t eend = time(NULL); + time_t fstart = time(NULL); + + if (mdtype == 1 || mdtype == 2) + { + mdb.runNVT(istep); + } + else if (mdtype == 0) + { + mdb.runNVE(istep); + } + else if (mdtype == -1) + { + stop = mdb.runFIRE(istep); + } + else + { + WARNING_QUIT("md_cells_classic", "mdtype should be -1~2!"); + } + time_t fend = time(NULL); + + ucell.save_cartesian_position(this->pos_next); + + ++istep; + } + + timer::tick("Run_MD_CLASSIC", "md_cells_classic"); return; +} + +void Run_MD_CLASSIC::md_allocate_ions(void) +{ + pos_dim = ucell.nat * 3; + + delete[] this->pos_old1; + delete[] this->pos_old2; + delete[] this->pos_now; + delete[] this->pos_next; + + this->pos_old1 = new double[pos_dim]; + this->pos_old2 = new double[pos_dim]; + this->pos_now = new double[pos_dim]; + this->pos_next = new double[pos_dim]; + + ZEROS(pos_old1, pos_dim); + ZEROS(pos_old2, pos_dim); + ZEROS(pos_now, pos_dim); + ZEROS(pos_next, pos_dim); +} + +void Run_MD_CLASSIC::update_pos_classic(void) +{ + for(int i=0;ipos_old2[i] = this->pos_old1[i]; + this->pos_old1[i] = this->pos_now[i]; + } + ucell.save_cartesian_position(this->pos_now); + return; } \ No newline at end of file diff --git a/source/module_md/run_md_classic.h b/source/module_md/run_md_classic.h index b1781d7dcd..841b102a98 100644 --- a/source/module_md/run_md_classic.h +++ b/source/module_md/run_md_classic.h @@ -2,6 +2,7 @@ #define RUN_MD_CLASSIC_H #include "../src_pw/tools.h" +#include "../module_neighbor/sltk_grid_driver.h" class Run_MD_CLASSIC { @@ -9,9 +10,20 @@ class Run_MD_CLASSIC Run_MD_CLASSIC(); ~Run_MD_CLASSIC(); + Grid_Driver grid_neigh; + void md_cells_classic(void); + void md_allocate_ions(void); + void update_pos_classic(void); + private: + int istep; + double* pos_old1; + double* pos_old2; + double* pos_now; + double* pos_next; + int pos_dim; }; diff --git a/source/run_lcao.cpp b/source/run_lcao.cpp index 51c433ce8d..60767242a2 100644 --- a/source/run_lcao.cpp +++ b/source/run_lcao.cpp @@ -8,6 +8,7 @@ #include "src_lcao/LOOP_cell.h" #include "src_io/print_info.h" #include "src_pw/symmetry.h" +#include "module_md/run_md_lcao.h" Run_lcao::Run_lcao(){} Run_lcao::~Run_lcao(){} @@ -143,8 +144,18 @@ void Run_lcao::lcao_line(void) } } - LOOP_cell lc; - lc.opt_cell(); + if(CALCULATION=="md") + { + Run_MD_LCAO run_md_lcao; + run_md_lcao.opt_cell(); + } + else // cell relaxations + { + LOOP_cell lc; + lc.opt_cell(); + + en.perform_dos(); + } en.perform_dos(); diff --git a/source/run_pw.cpp b/source/run_pw.cpp index 163331552d..5bfd99e7e5 100644 --- a/source/run_pw.cpp +++ b/source/run_pw.cpp @@ -9,6 +9,7 @@ #include "src_io/print_info.h" #include "src_pw/symmetry.h" #include "src_ions/Cell_PW.h" +#include "module_md/run_md_pw.h" Run_pw::Run_pw(){} Run_pw::~Run_pw(){} @@ -92,8 +93,16 @@ void Run_pw::plane_wave_line(void) CHR.allocate(NSPIN, pw.nrxx, pw.ngmc); pot.allocate(pw.nrxx); - Cell_PW cpws; - cpws.opt_cells_pw(); + if(CALCULATION == "md") + { + Run_MD_PW run_md_pw; + run_md_pw.md_cells_pw(); + } + else + { + Cell_PW cpws; + cpws.opt_cells_pw(); + } // caoyu add 2020-11-24, mohan updat 2021-01-03 diff --git a/source/src_io/print_info.cpp b/source/src_io/print_info.cpp index 76ac89a847..c4faf78bb2 100644 --- a/source/src_io/print_info.cpp +++ b/source/src_io/print_info.cpp @@ -107,11 +107,7 @@ void Print_Info::setup_parameters(void) cout << " ---------------------------------------------------------" << endl; - if(CALCULATION=="md" && INPUT.mdp.md_potential) - { - cout << " Classic Molecular Dynamics simulations" << endl; - } - else if(BASIS_TYPE=="lcao") + if(BASIS_TYPE=="lcao") { if(COLOUR && MY_RANK==0) { diff --git a/source/src_io/write_input.cpp b/source/src_io/write_input.cpp index 6f360d9eca..afe958443d 100644 --- a/source/src_io/write_input.cpp +++ b/source/src_io/write_input.cpp @@ -31,6 +31,7 @@ void Input::Print(const string &fn)const OUTP(ofs,"nbands_istate",nbands_istate,"number of bands around Fermi level for istate calulation"); OUTP(ofs,"nche_sto",nche_sto,"number of orders for Chebyshev expansion in stochastic DFT"); OUTP(ofs,"symmetry",symmetry,"turn symmetry on or off"); + OUTP(ofs,"set_vel",set_vel,"read velocity from STRU or not"); OUTP(ofs,"nelec",nelec,"input number of electrons"); ofs << "\n#Parameters (2.PW)" << endl; @@ -142,7 +143,7 @@ void Input::Print(const string &fn)const ofs << "\n#Parameters (10.Molecular dynamics)" << endl; OUTP(ofs,"md_mdtype",mdp.mdtype,"choose ensemble"); - OUTP(ofs,"md_potential",mdp.md_potential,"choose potential for md"); + //OUTP(ofs,"md_potential",mdp.md_potential,"choose potential for md"); OUTP(ofs,"md_dt",mdp.dt,"time step"); OUTP(ofs,"mnhc",mdp.MNHC,"number of Nose-Hoover chains"); OUTP(ofs,"md_qmass",mdp.Qmass,"mass of thermostat");