diff --git a/doc/generate-basis.md b/doc/generate-basis.md
index 29190dd75a0..33764bc3ab6 100644
--- a/doc/generate-basis.md
+++ b/doc/generate-basis.md
@@ -4,7 +4,7 @@
In ABACUS, the atomic orbital bases are generated using a scheme developed in the [paper](https://iopscience.iop.org/article/10.1088/0953-8984/22/44/445501). We provide a script named “generate_orbital.sh” under the directory tools/ to generate the atomic orbitals bases. In order to run this script, an ORBITAL_INPUT file is required.
-An example of this ORBITAL_INPUT file can be found in $ABACUS/tools/SIAB/2_Generate:
+An example of this ORBITAL_INPUT file can be found in $ABACUS/tools/SIAB/SimulatedAnnealing/example_N:
```
#1.exe_dir
#----------------------------------------------------------------------------
@@ -101,7 +101,7 @@ The ORBITAL_INPUT file contains 5 parts :
This part gives the bond lengths of the reference systems (dimer or trimer). Generally, the bond lengths are chosen to distribute on both sides of the equilibrium value. For example, for N dimer we use (in Å):
- Dis 1.0 1.1 1.5 2.0 3.0
- It means we take 5 reference systems (dimer), and the bond lengths are 1.0 1.1 1.5 2.0 3.0 angstrom, respectively. Every element has reference systems with different bond lengths, which could be found in file $ABACUS/tools/SIAB/2_Generate/DIS.txt.
+ It means we take 5 reference systems (dimer), and the bond lengths are 1.0 1.1 1.5 2.0 3.0 angstrom, respectively. Every element has reference systems with different bond lengths, which could be found in file $ABACUS/tools/SIAB/DIS.txt.
4. orbital generation
The main parameters for orbital generation
@@ -142,7 +142,7 @@ The ORBITAL_INPUT file contains 5 parts :
the accept rise of spillage when optimizing the kinetic energy
-After preparing the ORBITAL_INPUT file, one just needs to run the script and wait for the results. The results will be written into several output files under the directory $element.id_element/$Rcut/, for example 07_N/6/.
+After preparing the ORBITAL_INPUT file, one just needs to run the script "$PATH_TO/generate_orbital.sh ORBITAL_INPUT" and wait for the results. The results will be written into several output files under the directory $element.id_element/$Rcut/, for example 07_N/6/.
Some output files listed here are useful.
- ORBITAL_RESULTS.txt
@@ -160,4 +160,4 @@ For some elements, you can download the reference ORBITAL_INPUT files and pseudo
A file README is also given and you can decide the parameters with it as a reference.
In most cases, you just need to modify the parameters in Section 1, 2. Section 4 may be
partially modified if you need higher precision orbitals. The users are not encouraged to change
-the settings in sections 5, unless you are very familiar with the code generating algorithms.
\ No newline at end of file
+the settings in sections 5, unless you are very familiar with the code generating algorithms.
diff --git a/doc/input-main.md b/doc/input-main.md
index 4f9f06edc0d..c55542ef54f 100644
--- a/doc/input-main.md
+++ b/doc/input-main.md
@@ -19,15 +19,15 @@
- [Electronic structure](#electronic-structure)
- [basis_type](#basis-type) | [ks_solver](#ks-solver) | [nbands](#nbands) | [nbands_istate](#nbands-istate) | [nspin](#nspin) | [occupations](#occupations) | [smearing](#smearing) | [sigma](#sigma) | [mixing_type](#mixing-type) | [mixing_beta](#mixing-beta) | [mixing_ndim](#mixing-ndim) | [mixing_gg0](#mixing-gg0) | [gamma_only](#gamma-only) | [printe](#printe) | [niter](#niter) | [dr2](#dr2) | [charge_extrap](#charge-extrap)
+ [basis_type](#basis-type) | [ks_solver](#ks-solver) | [nbands](#nbands) | [nbands_istate](#nbands-istate) | [nspin](#nspin) | [occupations](#occupations) | [smearing](#smearing) | [sigma](#sigma) | [mixing_type](#mixing-type) | [mixing_beta](#mixing-beta) | [mixing_ndim](#mixing-ndim) | [mixing_gg0](#mixing-gg0) | [gamma_only](#gamma-only) | [printe](#printe) | [niter](#niter) | [dr2](#dr2) | [charge_extrap](#charge-extrap) | [ocp](#ocp) | [ocp_set](#ocp_set)
- [Geometry relaxation](#geometry-relaxation)
- [nstep](#nstep) | [force](#force) | [force_thr](#force-thr) | [force_thr_ev](#force-thr-ev) | [force_set](#force-set) | [bfgs_w1](#bfgs-w1) | [bfgs_w2](#bfgs-w2) | [trust_radius_max](#trust-radius-max) | [trust_radius_min](#trust-radius-min) | [trust_radius_ini](#trust-radius-ini) | [stress](#stress) | [stress_thr](#stress-thr) | [press](#press) | [fixed_axes](#fixed-axes) | [move_method](#move-method) | [cg_threshold](#cg-threshold) | [cell_factor](#cell-factor)
+ [nstep](#nstep) | [force](#force) | [force_thr](#force-thr) | [force_thr_ev](#force-thr-ev) | [force_set](#force-set) | [bfgs_w1](#bfgs-w1) | [bfgs_w2](#bfgs-w2) | [trust_radius_max](#trust-radius-max) | [trust_radius_min](#trust-radius-min) | [trust_radius_ini](#trust-radius-ini) | [stress](#stress) | [stress_thr](#stress-thr) | [press1, press2, press3](#press) | [fixed_axes](#fixed-axes) | [move_method](#move-method) | [cg_threshold](#cg-threshold) | [cell_factor](#cell-factor)
- [Variables related to program output](#variables-related-to-program-output)
- [mulliken](#mulliken) | [out_charge](#out-charge) | [out_potential](#out-potential) | [out_dm](#out-dm) | [out_wf](#out-wf) | [out_lowf](#out-lowf) | [out_dos](#out-dos) | [out_band](#out-band) | [out_stru](#out-stru) | [out_level](#out_level) | [out_alllog](#out-alllog) | [out_hs](#out-hs) | [out_r](#out-r) | [out_hs2](#out-hs2)
+ [mulliken](#mulliken) | [out_charge](#out-charge) | [out_potential](#out-potential) | [out_dm](#out-dm) | [out_wf](#out-wf) | [out_lowf](#out-lowf) | [out_dos](#out-dos) | [out_band](#out-band) | [out_stru](#out-stru) | [out_level](#out_level) | [out_alllog](#out-alllog) | [out_hs](#out-hs) | [out_r](#out-r) | [out_hs2](#out-hs2) | [out_element_info](#out-element-info) | [restart_save](#restart_save) | [restart_load](#restart_load)
- [Density of states](#density-of-states)
@@ -43,17 +43,23 @@
- [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)
+ [dft_plus_u](#dft_plus_u) | [orbital_corr](#orbital_corr) | [hubbard_u](#hubbard_u) | [hund_j](#hund_j) | [yukawa_potential](#yukawa_potential) | [omc](#omc)
+
- [VdW correction](#vdw-correction)
[vdw_method](#vdw-method) | [vdw_s6](#vdw-s6) | [vdw_s8](#vdw-s8) | [vdw_a1](#vdw-a1) | [vdw_a2](#vdw-a2) | [vdw_d](#vdw-d) | [vdw_abc](#vdw-abc) | [vdw_C6_file](#vdw-C6-file) | [vdw_C6_unit](#vdw-C6-unit) | [vdw_R0_file](#vdw-R0-file) | [vdw_R0_unit](#vdw-R0-unit) | [vdw_model](#vdw-model) | [vdw_radius](#vdw-radius) | [vdw_radius_unit](#vdw-radius-unit) | [vdw_cn_radius](#vdw-cn-radius) | [vdw_cn_radius_unit](#vdw-cn-radius-unit) | [vdw_period](#vdw-period)
- [Berry phase and wannier90 interface](#berry-phase-and-wannier90-interface)
- [berry_phase](#berry-phase) | [gdir](#gdir) | [towannier90](#towannier90) | [nnkpfile](#nnkpfile) | [wannier_spin](#wannier-spin) | [tddft](#tddft) [vext](#vext) | [vext_dire](#vext-dire)
+ [berry_phase](#berry-phase) | [gdir](#gdir) | [towannier90](#towannier90) | [nnkpfile](#nnkpfile) | [wannier_spin](#wannier-spin)
+
+ - [TDDFT: time dependent density functional theory](#TDDFT-doc)
+
+ [tddft](#tddft) | [td_dr2](#td_dr2) | [td_dt](#td_dt) | [td_force_dt](#td_force_dt) | [td_vext](#td_vext) | [td_vext_dire](#td_vext_dire) | [td_timescale](#td_timescale) | [td_vexttype](#td_vexttype) | [td_vextout](#td_vextout) | [td_dipoleout](#td_dipoleout)
- [Variables useful for debugging](#variables-useful-for-debugging)
@@ -292,7 +298,7 @@ This part of variables are used to control general system parameters.
- diago_proc
- *Type*: Integer
- - *Descrption*: If set to a positive number, then it specifies the number of threads used for carrying out diagonalization. Must be less than or equal to total number of MPI threads. Also, when cg diagonalization is used, diago_proc must be same as total number of MPI threads. If set to 0, then it will be set to the number of MPI threads. Normally, it is fine just leaving it to default value.
+ - *Descrption*: If set to a positive number, then it specifies the number of threads used for carrying out diagonalization. Must be less than or equal to total number of MPI threads. Also, when cg diagonalization is used, diago_proc must be same as total number of MPI threads. If set to 0, then it will be set to the number of MPI threads. Normally, it is fine just leaving it to default value. Only used for pw base.
- *Default*: 0
[back to top](#input-file)
@@ -605,6 +611,22 @@ calculations.
[back to top](#input-file)
+
+- ocp
+ - *Type*: Boolean
+ - *Description*: option for choose whether calcualting constrained DFT or not.
+ Only used for TDDFT.
+ - *Default*:0
+
+ [back to top](#input-file)
+
+- ocp_set
+ - *Type*: string
+ - *Description*: If ocp is true, the ocp_set is a string to set the number of occupancy, like 1 10 * 1 0 1 representing the 13 band occupancy, 12th band occupancy 0 and the rest 1, the code is parsing this string into an array through a regular expression.
+ - *Default*:none
+
+ [back to top](#input-file)
+
### Geometry relaxation
This part of variables are used to control the geometry relaxation.
@@ -687,11 +709,11 @@ This part of variables are used to control the geometry relaxation.
- stress_thr
- *Type*: Real
- *Description*: The threshold of the stress convergence, it indicates the largest stress among all the directions, the unit is KBar,
- - *Default*: 10
+ - *Default*: 0.01
[back to top](#input-file)
-- press1, 2, 3
+- press1, press2, press3
- *Type*: Real
- *Description*: the external pressures along three axes,the compressive stress is taken to be positive, the unit is KBar.
- *Default*: 0
@@ -831,6 +853,41 @@ This part of variables are used to control the output of properties.
[back to top](#input-file)
+- out_element_info
+ - *Type*: Boolean
+ - *Description*: When set to 1, ABACUS will generate a new directory under OUT.suffix path named as element name such as 'Si', which contained files "Si-d1-orbital-dru.dat Si-p2-orbital-k.dat Si-s2-orbital-dru.dat
+ Si-d1-orbital-k.dat Si-p2-orbital-r.dat Si-s2-orbital-k.dat
+ Si-d1-orbital-r.dat Si-p2-orbital-ru.dat Si-s2-orbital-r.dat
+ Si-d1-orbital-ru.dat Si-p-proj-k.dat Si-s2-orbital-ru.dat
+ Si.NONLOCAL Si-p-proj-r.dat Si-s-proj-k.dat
+ Si-p1-orbital-dru.dat Si-p-proj-ru.dat Si-s-proj-r.dat
+ Si-p1-orbital-k.dat Si-s1-orbital-dru.dat Si-s-proj-ru.dat
+ Si-p1-orbital-r.dat Si-s1-orbital-k.dat v_loc_g.dat
+ Si-p1-orbital-ru.dat Si-s1-orbital-r.dat
+ Si-p2-orbital-dru.dat Si-s1-orbital-ru.dat " for example.
+
+ - *Default*: 0
+
+ [back to top](#input-file)
+
+- restart_save
+ - *Type*: Boolean
+ - *Description*: Only for LCAO, store charge density file and H matrix file every scf step for restart.
+ - *Default*: 0
+
+ [back to top](#input-file)
+
+- restart_load
+ - *Type*: Boolean
+ - *Description*: Only for LCAO, used for restart, only if that:
+ * set restart_save as true and do scf calculation before.
+ * please ensure suffix is same with calculation before and density file and H matrix file is exist.
+
+ restart from stored density file and H matrix file.
+ - *Default*: 0
+
+ [back to top](#input-file)
+
### Density of states
This part of variables are used to control the calculation of DOS.
@@ -859,7 +916,7 @@ This part of variables are used to control the calculation of DOS.
This part of variables are used to control the addition of an external electric field. It is achieved by adding a saw-like potential to the local ionic potential.
- efield
- - *Type*: Bool
+ - *Type*: Boolean
- *Description*: Controls whether to add the external electric field. When set to 1, the electric field is turned on. When set to 0, there is no electric field.
- *Default*: 0.
@@ -895,9 +952,10 @@ This part of variables are used to control the addition of an external electric
### DeePKS
This part of variables are used to control the usage of DeePKS method (a comprehensive data-driven approach to improve accuracy of DFT).
+Warning: this function is not robust enough for version 2.2.0. Please try these variables in https://github.com/deepmodeling/abacus-develop/tree/deepks .
- out_descriptor
- - *Type*: Bool
+ - *Type*: Boolean
- *Description*: when set to 1, ABACUS will calculate and output descriptor for DeePKS training. In `LCAO` calculation, a path of *.orb file is needed to be specified under `NUMERICAL_DESCRIPTOR`in `STRU`file. For example:
```
NUMERICAL_ORBITAL
@@ -917,7 +975,7 @@ This part of variables are used to control the usage of DeePKS method (a compreh
[back to top](#input-file)
- deepks_scf
- - *Type*: Bool
+ - *Type*: Boolean
- *Description*: only when deepks is enabled in `LCAO` calculation can this variable set to 1. Then, a trained, traced model file is needed for self-consistant field iteration in DeePKS method.
- *Default*: 0
@@ -1053,10 +1111,13 @@ This part of variables are used to control the molecular dynamics calculations.
- md_type
- *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)
@@ -1071,17 +1132,17 @@ This part of variables are used to control the molecular dynamics calculations.
[back to top](#input-file)
- md_rstmd
- - *Type*: Bool
+ - *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
- *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
@@ -1092,8 +1153,8 @@ This part of variables are used to control the molecular dynamics calculations.
[back to top](#input-file)
- md_qmass
- *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
@@ -1103,6 +1164,14 @@ This part of variables are used to control the molecular dynamics calculations.
[back to top](#input-file)
+- md_tfreq
+ - *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
- *Type*: Integer
@@ -1116,7 +1185,7 @@ This part of variables are used to control the molecular dynamics calculations.
- NVT_control
- *Type*: Integer
- *Description*: Specifies which type of thermostat is used.
- - 1: Nose-Hoover
+ - 1: Nose-Hoover-chains
- 2: Langevin
- 3: Andersen
- *Default*: 1
@@ -1172,45 +1241,80 @@ This part of variables are used to control the molecular dynamics calculations.
[back to top](#input-file)
+- direction
+ - *Type*: Integer
+ - *Description*: the direction of shock wave for MSST.
+ - *Default*: 2 (z direction)
+
+ [back to top](#input-file)
+
+- velocity
+ - *Type*: Real
+ - *Description*: the velocity of shock wave (\AA/fs) for MSST.
+ - *Default*: 0
+
+ [back to top](#input-file)
+
+- viscosity
+ - *Type*: Real
+ - *Description*: artificial viscosity (mass/length/time) for MSST.
+ - *Default*: 0
+
+ [back to top](#input-file)
+
+- tscale
+ - *Type*: Real
+ - *Description*: reduction in initial temperature (0~1) used to compress volume in MSST.
+ - *Default*: 0
+
+ [back to top](#input-file)
+
+- damp
+ - *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
- - *Type*: Bool
+- dft_plus_u
+ - *Type*: Boolean
- *Description*: If set to 1, ABCUS will calculate plus U correction, which is especially important for correlated electron.
- *Default*: 0
[back to top](#input-file)
-- orbital_corr
+- orbital_corr
- *Type*: Int
- *Description*: $l_1,l_2,l_3,\ldots$ for atom type 1,2,3 respectively.(usually 2 for d electrons and 3 for f electrons) .Specify which orbits need plus U correction for each atom. If set to -1, the correction would not be calculate for this atom.
- *Default*: None
[back to top](#input-file)
-- hubbard_u
+- hubbard_u
- *Type*: Real
- *Description*: Hubbard Coulomb interaction parameter U(ev) in plus U correction,which should be specified for each atom unless Yukawa potential is use. ABACUS use a simplified scheme which only need U and J for each atom.
- *Default*: 0.0
[back to top](#input-file)
-- hund_j
+- hund_j
- *Type*: Real
- *Description*: Hund exchange parameter J(ev) in plus U correction ,which should be specified for each atom unless Yukawa potential is use. ABACUS use a simplified scheme which only need U and J for each atom.
- *Default*: 0.0
[back to top](#input-file)
-- yukawa_potential
- - *Type*: Bool
+- yukawa_potential
+ - *Type*: Boolean
- *Description*: whether use the local screen Coulomb potential method to calculate the value of U and J. If this is set to 1, hubbard_u and hund_j do not need to be specified.
- *Default*: 0
[back to top](#input-file)
-- omc
- - *Type*: Bool
+- omc
+ - *Type*: Boolean
- *Description*: whether turn on occupation matrix control method or not
- *Default*: 0
@@ -1361,6 +1465,8 @@ This part of variables are used to control berry phase and wannier90 interfacae
- *Default*: up
[back to top](#input-file)
+
+### TDDFT: time dependent density functional theory
- tddft
- *Type*: Integer
- *Description*:
@@ -1369,7 +1475,25 @@ This part of variables are used to control berry phase and wannier90 interfacae
- *Default*: 0
[back to top](#input-file)
-- vext
+- td_dr2
+ - *Type*: Double
+ - *Description*: Accuracy of electron convergence when doing time-dependent evolution.
+ - *Default*: 1e-9
+
+ [back to top](#input-file)
+- td_dt
+ - *Type*: Double
+ - *Description*: Time-dependent evolution time step. (fs)
+ - *Default*: 0.02
+
+ [back to top](#input-file)
+- td_force_dt
+ - *Type*: Double
+ - *Description*: Time-dependent evolution force changes time step. (fs)
+ - *Default*: 0.02
+
+ [back to top](#input-file)
+- td_vext
- *Type*: Integer
- *Description*:
- 1: add a laser material interaction (extern laser field).
@@ -1377,7 +1501,7 @@ This part of variables are used to control berry phase and wannier90 interfacae
- *Default*: 0
[back to top](#input-file)
-- vext_dire
+- td_vext_dire
- *Type*: Integer
- *Description*:
- 1: the direction of external light field is along x axis.
@@ -1385,6 +1509,37 @@ This part of variables are used to control berry phase and wannier90 interfacae
- 3: the direction of external light field is along z axis.
- *Default*: 1
+ [back to top](#input-file)
+- td_timescale
+ - *Type*: Double
+ - *Description*: Time range of external electric field application. (fs)
+ - *Default*: 0.5
+
+ [back to top](#input-file)
+- td_vexttype
+ - *Type*: Integer
+ - *Description*:
+ - 1: Gaussian-type light field.
+ - 2: Delta function form light field.
+ - 3: Trigonometric function form light field.
+ - *Default*: 1
+
+ [back to top](#input-file)
+- td_vextout
+ - *Type*: Integer
+ - *Description*:
+ - 1: Output external electric field.
+ - 0: do not Output external electric field.
+ - *Default*: 0
+
+ [back to top](#input-file)
+- td_dipoleout
+ - *Type*: Integer
+ - *Description*:
+ - 1: Output dipole.
+ - 0: do not Output dipole.
+ - *Default*: 0
+
[back to top](#input-file)
### Variables useful for debugging
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index a7ff0bc5f66..4df01cc1cb6 100644
--- a/source/Makefile.Objects
+++ b/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\
@@ -256,6 +253,7 @@ write_rho_dipole.o\
write_HS.o\
write_HS_R.o\
write_dm.o\
+write_wfc_realspace.o\
potential_libxc.o \
potential_libxc_meta.o \
efield.o \
@@ -281,6 +279,7 @@ variable_cell.o\
dftu.o\
dftu_yukawa.o\
dftu_relax.o\
+dmft.o \
OBJS_COMMON=atom_spec.o \
unitcell.o \
@@ -297,6 +296,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 \
diff --git a/source/input.cpp b/source/input.cpp
index 773a49b6078..a87ad57ddd1 100644
--- a/source/input.cpp
+++ b/source/input.cpp
@@ -137,7 +137,6 @@ void Input::Default(void)
opt_epsilon2 = false;//mohan add 2010-03-24
opt_nbands = 0;
- lda_plus_u = false;
//----------------------------------------------------------
// electrons / spin
//----------------------------------------------------------
@@ -158,7 +157,6 @@ void Input::Default(void)
symmetry=false;
set_vel=false;
symmetry_prec = 1.0e-5; //LiuXh add 2021-08-12, accuracy for symmetry
- mlwf_flag=false;
force=0;
force_set=false;
force_thr=1.0e-3;
@@ -221,7 +219,7 @@ void Input::Default(void)
//----------------------------------------------------------
dr2 = 1.0e-9;
niter = 40;
- nstep = 1;
+ this->nstep = 0;
out_stru = 0;
//----------------------------------------------------------
// occupation
@@ -253,6 +251,7 @@ void Input::Default(void)
out_potential = 0;
out_wf = 0;
+ out_wf_r = 0;
out_dos = 0;
out_band = 0;
out_hs = 0;
@@ -265,6 +264,7 @@ void Input::Default(void)
dos_edelta_ev = 0.01;//(ev)
dos_scale = 0.01;
b_coef = 0.07;
+ out_element_info = false;
//----------------------------------------------------------
// LCAO
//----------------------------------------------------------
@@ -282,33 +282,6 @@ void Input::Default(void)
selinv_mu = -1.0;
selinv_threshold = 1.0e-3;
selinv_niter = 50;
-//----------------------------------------------------------
-// Molecular Dynamics
-//----------------------------------------------------------
-/*
- md_dt=20.0; //unit is 1 a.u., which is 4.8378*10e-17 s
- md_restart=0;
- md_tolv=100.0;
- md_thermostat="not_controlled"; //"rescaling","rescale-v","rescale-t","reduce-t"...
- md_temp0=300; //kelvin
- md_tstep=1; //reduec md_delt every md_tstep step.
- md_delt=1.0;
-*/
-
-/* //----------------------------------------------------------
-// vdwD2 //Peize Lin add 2014-03-31, update 2015-09-30
-//----------------------------------------------------------
- vdwD2=false;
- vdwD2_scaling=0.75;
- vdwD2_d=20;
- vdwD2_C6_file="default";
- vdwD2_C6_unit="Jnm6/mol";
- vdwD2_R0_file="default";
- vdwD2_R0_unit="A";
- vdwD2_model="radius";
- vdwD2_period = {3,3,3};
- vdwD2_radius=30.0/ModuleBase::BOHR_TO_A;
- vdwD2_radius_unit="Bohr"; */
//----------------------------------------------------------
// vdw //jiyy add 2019-08-04
@@ -466,6 +439,11 @@ void Input::Default(void)
omc = false;
dftu_type = 2;
+//==========================================================
+// DFT+DMFT Xin Qu added on 2020-08
+//==========================================================
+ dft_plus_dmft = false;
+
return;
}
@@ -663,10 +641,6 @@ bool Input::Read(const std::string &fn)
{
read_value(ifs, opt_nbands);
}
- else if (strcmp("lda_plus_u", word) == 0)// lda + u
- {
- read_value(ifs, lda_plus_u);
- }
//----------------------------------------------------------
// electrons / spin
//----------------------------------------------------------
@@ -730,10 +704,6 @@ bool Input::Read(const std::string &fn)
{
read_value(ifs, symmetry_prec);
}
- else if (strcmp("mlwf_flag", word) == 0)
- {
- read_value(ifs, mlwf_flag);
- }
else if (strcmp("force", word) == 0)
{
read_value(ifs, force);
@@ -947,7 +917,7 @@ bool Input::Read(const std::string &fn)
}
else if (strcmp("nstep", word) == 0)
{
- read_value(ifs, nstep);
+ read_value(ifs, this->nstep);
}
else if (strcmp("out_stru", word) == 0)
{
@@ -1060,6 +1030,10 @@ bool Input::Read(const std::string &fn)
else if (strcmp("out_wf", word) == 0)
{
read_value(ifs, out_wf);
+ }
+ else if (strcmp("out_wf_r", word) == 0)
+ {
+ read_value(ifs, out_wf_r);
}
//mohan add 20090909
else if (strcmp("out_dos", word) == 0)
@@ -1091,6 +1065,10 @@ bool Input::Read(const std::string &fn)
else if (strcmp("out_alllog", word) == 0)
{
read_value(ifs, out_alllog);
+ }
+ else if (strcmp("out_element_info", word) == 0)
+ {
+ read_value(ifs, out_element_info);
}
else if (strcmp("dos_emin_ev", word) == 0)
{
@@ -1162,37 +1140,7 @@ bool Input::Read(const std::string &fn)
read_value(ifs, selinv_niter);
}
// about molecular dynamics
-/*
- else if (strcmp("md_dt", word) == 0)
- {
- read_value(ifs, md_dt);
- }
- else if (strcmp("md_restart", word) == 0)
- {
- read_value(ifs, md_restart);
- }
- else if (strcmp("md_tolv", word) == 0)
- {
- read_value(ifs, md_tolv);
- }
- else if (strcmp("md_thermostat", word) == 0)
- {
- read_value(ifs, md_thermostat);
- }
- else if (strcmp("md_temp0", word) == 0)
- {
- read_value(ifs, md_temp0);
- }
- else if (strcmp("md_tstep", word) == 0)
- {
- read_value(ifs, md_tstep);
- }
- else if (strcmp("md_delt", word) == 0)
- {
- read_value(ifs, md_delt);
- }
-*/
-//added begin by zheng daye
+ //added begin by zheng daye
else if (strcmp("md_mdtype",word) == 0)
{
read_value(ifs, mdp.mdtype);
@@ -1249,7 +1197,7 @@ bool Input::Read(const std::string &fn)
{
read_value(ifs,mdp.ediffg );
}
-//added by zheng daye
+ //added by zheng daye
//----------------------------------------------------------
// Classic MD
// Yu Liu add 2021-07-30
@@ -1270,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
@@ -1326,57 +1298,6 @@ bool Input::Read(const std::string &fn)
{
read_value(ifs,td_dipoleout );
}
-
-
-/* //----------------------------------------------------------
-// vdwD2
-// Peize Lin add 2014-03-31
-//----------------------------------------------------------
- else if (strcmp("vdwd2", word) == 0)
- {
- read_value(ifs, vdwD2);
- }
- else if (strcmp("vdwd2_scaling", word) == 0)
- {
- read_value(ifs, vdwD2_scaling);
- }
- else if (strcmp("vdwd2_d", word) == 0)
- {
- read_value(ifs, vdwD2_d);
- }
- else if (strcmp("vdwd2_c6_file", word) == 0)
- {
- read_value(ifs, vdwD2_C6_file);
- }
- else if (strcmp("vdwd2_c6_unit", word) == 0)
- {
- read_value(ifs, vdwD2_C6_unit);
- }
- else if (strcmp("vdwd2_r0_file", word) == 0)
- {
- read_value(ifs, vdwD2_R0_file);
- }
- else if (strcmp("vdwd2_r0_unit", word) == 0)
- {
- read_value(ifs, vdwD2_R0_unit);
- }
- else if (strcmp("vdwd2_model", word) == 0)
- {
- read_value(ifs, vdwD2_model);
- }
- else if (strcmp("vdwd2_period", word) == 0)
- {
- ifs >> vdwD2_period.x >> vdwD2_period.y;
- read_value(ifs, vdwD2_period.z);
- }
- else if (strcmp("vdwd2_radius", word) == 0)
- {
- read_value(ifs, vdwD2_radius);
- }
- else if (strcmp("vdwd2_radius_unit", word) == 0)
- {
- read_value(ifs, vdwD2_radius_unit);
- } */
//----------------------------------------------------------
// vdw
// jiyy add 2019-08-04
@@ -1464,14 +1385,6 @@ bool Input::Read(const std::string &fn)
//--------------------------------------------------------
// epsilon pengfei Li 2016-11-23
//--------------------------------------------------------
- //else if (strcmp("epsilon", word) == 0)
- //{
- // read_value(ifs, epsilon);
- //}
- //else if (strcmp("epsilon_choice", word) == 0)
- //{
- // read_value(ifs, epsilon_choice);
- //}
else if (strcmp("spectral_type", word) == 0)
{
read_value(ifs, spectral_type);
@@ -1512,10 +1425,6 @@ bool Input::Read(const std::string &fn)
{
read_value(ifs, ecut_chi);
}
- //else if (strcmp("oband", word) == 0)
- //{
- // read_value(ifs, oband);
- //}
else if (strcmp("q_start", word) == 0)
{
ifs >> q_start[0]; ifs >> q_start[1]; read_value(ifs, q_start[2]);
@@ -1524,14 +1433,6 @@ bool Input::Read(const std::string &fn)
{
ifs >> q_direct[0]; ifs >> q_direct[1]; read_value(ifs, q_direct[2]);
}
- //else if (strcmp("start_q", word) == 0)
- //{
- // read_value(ifs, start_q);
- //}
- //else if (strcmp("interval_q", word) == 0)
- //{
- // read_value(ifs, interval_q);
- //}
else if (strcmp("nq", word) == 0)
{
read_value(ifs, nq);
@@ -1576,18 +1477,6 @@ bool Input::Read(const std::string &fn)
getline(ifs, GlobalV::ocp_set);
// ifs.ignore(150, '\n');
}
- // else if (strcmp("ocp_n", word) == 0)
- // {
- // read_value(ifs, ocp_n);
- // }
- // else if (strcmp("ocp_kb", word) == 0)
- // {
- // for(int i=0; i<(ocp_n-1); i++)
- // {
- // ifs >> GlobalV::ocp_kb[i];
- // }
- // read_value(ifs, GlobalV::ocp_kb[ocp_n-1]);
- // }
else if (strcmp("mulliken", word) == 0)
{
read_value(ifs, GlobalV::mulliken);
@@ -1597,14 +1486,6 @@ bool Input::Read(const std::string &fn)
ifs >> lcao_box[0]; ifs >> lcao_box[1];
read_value(ifs, lcao_box[2]);
}
- //else if (strcmp("epsilon0", word) == 0)
- //{
- // read_value(ifs, epsilon0);
- //}
- //else if (strcmp("intersmear", word) == 0)
- //{
- // read_value(ifs, intersmear);
- //}
else if (strcmp("intrasmear", word) == 0)
{
read_value(ifs, intrasmear);
@@ -1621,10 +1502,6 @@ bool Input::Read(const std::string &fn)
{
read_value(ifs, eps_degauss);
}
- //else if (strcmp("epsilon0_choice", word) == 0)
- //{
- // read_value(ifs, epsilon0_choice);
- //}
//----------------------------------------------------------
// exx
// Peize Lin add 2018-06-20
@@ -1713,22 +1590,6 @@ bool Input::Read(const std::string &fn)
{
read_value(ifs, soc_lambda);
}
-/* else if (strcmp("angle1", word) == 0)
- {
- angle1.resize(ntype);
- for(auto &i:angle1)
- read_value(ifs, i);
- }
- else if (strcmp("angle2", word) == 0)
- {
- angle2.resize(ntype);
- for (auto &i : angle2)
- read_value(ifs, i);
- }*/
- //else if (strcmp("epsilon0_choice", word) == 0)
- //{
- // read_value(ifs, epsilon0_choice);
- //}
else if (strcmp("cell_factor", word) == 0)
{
read_value(ifs, cell_factor);
@@ -1741,81 +1602,6 @@ bool Input::Read(const std::string &fn)
{
read_value(ifs, test_just_neighbor);
}
-//---------------
-//start magnetic
-/*
-#ifndef __CMD
- else if (strcmp("magmom", word) == 0)
- {
- n_mag_at=0;
- stringstream sstr;
- string s;
- getline(ifs,s);
- sstr.str(s);
- int tmplength=s.length();
- int at_per_mag[tmplength];
- double mags[tmplength];
- int n_magmom=0;
- s="";
- //1 3*2 6*1
- while(sstr.good())
- {
- sstr>>s;
- string s1;
- string s2;
- bool mul=0;// if this parameter is the form n*m
- for(int i=0;i='0'&& s[i]<='9') or s[i]=='.' or s[i]=='+' or s[i]=='-')
- {
- s1.push_back(s[i]);
- }
- else if (s[i]=='*')
- {
- s2=s1;
- s1="";
- mul=true;
- }
- else
- {
- std::cout<<"Unrecognized character"<=at_per_mag[n_m])
- {
- n_n+=at_per_mag[n_m];
- n_m+=1;
- }
- atom_mag[i]=mags[n_m];
- cout<<"atom_mag"<> dft_plus_dmft;
+ }
+//----------------------------------------------------------------------------------
else
{
//xiaohui add 2015-09-15
@@ -1978,7 +1771,7 @@ bool Input::Read(const std::string &fn)
exit(0);
}
- if( (orbital_corr[i]==-1) && (orbital_corr[i]==0) && (orbital_corr[i]!=1) && (orbital_corr[i]!=2) && (orbital_corr[i]!=3) )
+ if( (orbital_corr[i]!=-1) && (orbital_corr[i]!=0) && (orbital_corr[i]!=1) && (orbital_corr[i]!=2) && (orbital_corr[i]!=3) )
{
std::cout << " WRONG ARGUMENTS OF orbital_corr " << std::endl;
exit(0);
@@ -1997,9 +1790,102 @@ bool Input::Read(const std::string &fn)
exit(0);
}
+ if(strcmp("genelpa", ks_solver.c_str())!=0 && strcmp(ks_solver.c_str(),"scalapack_gvx")!=0 )
+ {
+ std::cout << " WRONG ARGUMENTS OF ks_solver in DFT+U routine, only genelpa and scalapack_gvx are supportted " << std::endl;
+ exit(0);
+ }
+
+ }
+
+//----------------------------------------------------------
+// DFT+DMFT Xin Qu added on 2020-08
+//----------------------------------------------------------
+ if(dft_plus_dmft)
+ {
+ ifs.clear();
+ ifs.seekg(0); //move to the beginning of the file
+ ifs.rdstate();
+ while (ifs.good())
+ {
+ ifs >> word1;
+ strtolower(word1, word); //convert uppercase std::string to lower case; word1 --> word
+
+ if(strcmp("hubbard_u", word)==0)
+ {
+ for(int i=0; i> hubbard_u[i];
+ hubbard_u[i] /= ModuleBase::Ry_to_eV;
+ }
+ }
+ else if (strcmp("hund_j", word)==0)
+ {
+ for(int i=0;i> hund_j[i];
+ hund_j[i] /= ModuleBase::Ry_to_eV;
+ }
+ }
+ else if(strcmp("orbital_corr", word)==0)
+ {
+ for(int i=0;i> orbital_corr[i];
+ }
+ }
+ else ifs.ignore(150, '\n');
+
+ if (ifs.eof() != 0) break;
+ }
+
+ for(int i=0; instep );
Parallel_Common::bcast_int( out_stru ); //mohan add 2012-03-23
//Parallel_Common::bcast_string( occupations );
@@ -2227,6 +2111,7 @@ void Input::Bcast()
Parallel_Common::bcast_int(out_potential);
Parallel_Common::bcast_int( out_wf );
+ Parallel_Common::bcast_int( out_wf_r );
Parallel_Common::bcast_int( out_dos );
Parallel_Common::bcast_int( out_band );
Parallel_Common::bcast_int( out_hs );
@@ -2234,6 +2119,7 @@ void Input::Bcast()
Parallel_Common::bcast_int( out_r_matrix ); // jingan add 2019-8-14
Parallel_Common::bcast_bool( out_lowf );
Parallel_Common::bcast_bool( out_alllog );
+ Parallel_Common::bcast_bool( out_element_info );
Parallel_Common::bcast_double( dos_emin_ev );
Parallel_Common::bcast_double( dos_emax_ev );
@@ -2284,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 );
@@ -2439,7 +2331,7 @@ void Input::Bcast()
//-----------------------------------------------------------------------------------
//DFT+U (added by Quxin 2020-10-29)
//-----------------------------------------------------------------------------------
- Parallel_Common::bcast_bool( dft_plus_u );
+ Parallel_Common::bcast_bool( dft_plus_u );
Parallel_Common::bcast_bool( yukawa_potential );
Parallel_Common::bcast_bool( omc );
Parallel_Common::bcast_int(dftu_type);
@@ -2459,6 +2351,11 @@ void Input::Bcast()
Parallel_Common::bcast_int(orbital_corr[i]);
}
+//-----------------------------------------------------------------------------------
+//DFT+DMFT (added by Quxin 2020-08)
+//-----------------------------------------------------------------------------------
+ Parallel_Common::bcast_bool( dft_plus_dmft );
+
return;
}
#endif
@@ -2475,6 +2372,10 @@ void Input::Check(void)
//std::cout << "diago_proc=" << diago_proc << std::endl;
//std::cout << " NPROC=" << GlobalV::NPROC << std::endl;
+ if(diago_proc>1 && basis_type=="lcao")
+ {
+ ModuleBase::WARNING_QUIT("Input", "please don't set diago_proc with lcao base");
+ }
if(diago_proc<=0)
{
diago_proc = GlobalV::NPROC;
@@ -2484,19 +2385,6 @@ void Input::Check(void)
diago_proc = GlobalV::NPROC;
}
- // mohan add 2010/03/29
- //if(!local_basis && diago_type=="lapack") xiaohui modify 2013-09-01
- //if(basis_type=="pw" && ks_solver=="lapack") xiaohui modify 2013-09-04 //xiaohui add 2013-09-01
- //{
- // ModuleBase::WARNING_QUIT("Input","lapack can not be used in plane wave basis.");
- //} xiaohui modify 2013-09-04
-
- //xiaohui move 4 lines, 2015-09-30
- //if(symmetry)
- //{
- // ModuleBase::WARNING("Input","symmetry is only correct for total energy calculations now,not for nonlocal force." );
- //}
-
if (efield && symmetry)
{
symmetry = false;
@@ -2556,7 +2444,7 @@ void Input::Check(void)
std::cout<<"sorry, can't calculate force with soc now, would be implement in next version!"<nstep = 1;
}
else if (calculation == "scf-sto") // qianrui 2021-2-20
@@ -2566,6 +2454,7 @@ void Input::Check(void)
mem_saver = 0;
ModuleBase::GlobalFunc::AUTO_SET("mem_savre","0");
}
+ this->nstep = 1;
}
else if (calculation == "relax") // pengfei 2014-10-13
{
@@ -2575,12 +2464,13 @@ void Input::Check(void)
ModuleBase::GlobalFunc::AUTO_SET("mem_savre","0");
}
force = 1;
+ if(! this->nstep) this->nstep = 50;
}
else if (calculation == "nscf")
{
GlobalV::CALCULATION = "nscf";
- nstep = 1;
+ this->nstep = 1;
out_stru = 0;
//if (local_basis == 0 && linear_scaling == 0) xiaohui modify 2013-09-01
@@ -2604,7 +2494,7 @@ void Input::Check(void)
else if(calculation == "istate")
{
GlobalV::CALCULATION = "istate";
- nstep = 1;
+ this->nstep = 1;
out_stru = 0;
out_dos = 0;
out_band = 0;
@@ -2625,7 +2515,7 @@ void Input::Check(void)
else if(calculation == "ienvelope")
{
GlobalV::CALCULATION = "ienvelope"; // mohan fix 2011-11-04
- nstep = 1;
+ this->nstep = 1;
out_stru = 0;
out_dos = 0;
out_band = 0;
@@ -2647,6 +2537,10 @@ void Input::Check(void)
GlobalV::CALCULATION = "md";
symmetry = false;
force = 1;
+ if(this->nstep==0){
+ GlobalV::ofs_running<<"nstep should be set. Autoset nstep to 50!"<nstep = 50;
+ }
if(!out_md_control) out_level = "m";//zhengdy add 2019-04-07
//deal with input parameters , 2019-04-30
@@ -2677,9 +2571,11 @@ void Input::Check(void)
{
force = 1;
stress = 1;
+ if(! this->nstep) this->nstep = 50;
}
else if(calculation == "test")
{
+ this->nstep = 1;
}
else
{
@@ -3005,14 +2901,6 @@ void Input::Check(void)
{
ModuleBase::WARNING_QUIT("INPUT","You must choose a direction!");
}
- //if( oband > nbands)
- //{
- // ModuleBase::WARNING_QUIT("INPUT","oband must <= nbands");
- //}
- // if( oband == 1)
- // {
- // oband = nbands;
- // }
}
if(exx_hybrid_type!="no" &&
diff --git a/source/input.h b/source/input.h
index 5d1c3eab523..7c8047435e6 100644
--- a/source/input.h
+++ b/source/input.h
@@ -63,8 +63,6 @@ class Input
std::string NNKP; // add by jingan for wannier90
std::string wannier_spin; // add by jingan for wannier90
- bool mlwf_flag; // add by mohan
-
//==========================================================
// Stochastic DFT
//==========================================================
@@ -89,7 +87,6 @@ class Input
//==========================================================
bool opt_epsilon2; // true : calculate the dielectric functions
int opt_nbands; // number of bands for optical transition matrix
- bool lda_plus_u; // true : lda plus u calculation
//==========================================================
// electrons / spin
@@ -222,6 +219,7 @@ class Input
int out_dm; // output density matrix.
int out_potential; // yes or no
int out_wf; // 0: no; 1: txt; 2: dat
+ int out_wf_r; // 0: no; 1: yes
int out_dos; // dos calculation. mohan add 20090909
int out_band; // band calculation pengfei 2014-10-13
int out_hs; // output H matrix and S matrix in local basis.
@@ -229,6 +227,7 @@ class Input
int out_r_matrix; // jingan add 2019-8-14, output r(R) matrix.
bool out_lowf; // output the wave functions in local basis.
bool out_alllog; // output all logs.
+ bool out_element_info; // output infomation of all element
double dos_emin_ev;
double dos_emax_ev;
@@ -447,6 +446,11 @@ class Input
int dftu_type; //1:rotationally invarient formalism; 2:simplified form(default)
int double_counting; // 1:FLL(fully localized limit)(default); 2:AMF(around mean field)
+//==========================================================
+// DFT+DMFT Xin Qu added on 2021-08
+//==========================================================
+ bool dft_plus_dmft; //true:DFT+U correction; false:standard DFT calcullation(default)
+
//==========================================================
// DeepKS -- added by caoyu and mohan
//==========================================================
diff --git a/source/input_conv.cpp b/source/input_conv.cpp
index 1755ebd37d6..5f9f1d1f532 100644
--- a/source/input_conv.cpp
+++ b/source/input_conv.cpp
@@ -83,6 +83,7 @@ void Input_Conv::Convert(void)
GlobalV::PRESS1 = INPUT.press1;
GlobalV::PRESS2 = INPUT.press2;
GlobalV::PRESS3 = INPUT.press3;
+ GlobalV::out_element_info = INPUT.out_element_info;
#ifdef __LCAO
Force_Stress_LCAO::force_invalid_threshold_ev = INPUT.force_thr_ev2;
#endif
@@ -398,8 +399,10 @@ void Input_Conv::Convert(void)
- // jiyy add 2020.10.11
- // fix bugs of ocp_set -- Yuanbo Li 2021/8/17
+ // setting for constrained DFT, jiyy add 2020.10.11
+ // For example, when we studying nitrogen-vacancy center,
+ // it requires an additional excitation of an electron conduction band to simulate the excited state,
+ // used for TDDFT only.
if(GlobalV::ocp == 1)
{
int count = 0;
@@ -591,6 +594,7 @@ void Input_Conv::Convert(void)
GlobalC::CHR.nelec = INPUT.nelec;
GlobalC::pot.out_potential = INPUT.out_potential;
GlobalC::wf.out_wf = INPUT.out_wf;
+ GlobalC::wf.out_wf_r = INPUT.out_wf_r;
GlobalC::en.out_dos = INPUT.out_dos;
GlobalC::en.out_band = INPUT.out_band;
#ifdef __LCAO
diff --git a/source/module_base/constants.h b/source/module_base/constants.h
index d863028e561..86e10f44e6b 100644
--- a/source/module_base/constants.h
+++ b/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
diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp
index 31983fadbf4..f1bb46a1604 100644
--- a/source/module_base/global_variable.cpp
+++ b/source/module_base/global_variable.cpp
@@ -185,4 +185,6 @@ bool out_descriptor = false; //caoyu add 2021-10-16 for DeePKS
bool deepks_scf = false; //caoyu add 2021-10-16 for DeePKS
int vnl_method = 1; //set defauld vnl method as old, added by zhengdy 2021-10-11
+
+bool out_element_info = false; //added by zhengdy 2021-11-26
}
diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h
index 75f9373d4e5..9c68ae56893 100644
--- a/source/module_base/global_variable.h
+++ b/source/module_base/global_variable.h
@@ -203,5 +203,7 @@ extern bool deepks_scf; //caoyu add 2021-10-16 for DeePKS
//method for dealing with non-local potential in Hamiltonian matrix, 0 for old, 1 for new
extern int vnl_method;
+//whether or not output information for each element
+extern bool out_element_info;
}
#endif
diff --git a/source/module_base/matrix3.h b/source/module_base/matrix3.h
index 9167e99be16..5a120b5926b 100644
--- a/source/module_base/matrix3.h
+++ b/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 */
diff --git a/source/module_base/timer.cpp b/source/module_base/timer.cpp
index 22ac495546c..cf29c267be4 100644
--- a/source/module_base/timer.cpp
+++ b/source/module_base/timer.cpp
@@ -149,12 +149,13 @@ void timer::print_all(std::ofstream &ofs)
if(timer_one.cpu_second < small)
continue;
+ ofs << std::resetiosflags(std::ios::scientific);
ofs << " "
// << std::setw(2) << timer_one.level
<< std::setw(2) << " "
<< std::setw(20) << class_name
<< std::setw(20) << name
- << std::setw(15) << timer_one.cpu_second
+ << std::setw(15) << std::setprecision(5) << timer_one.cpu_second
<< std::setw(10) << timer_one.calls
<< std::setw(10) << std::setprecision(2) << timer_one.cpu_second/timer_one.calls
<< std::setw(10) << timer_one.cpu_second / timer_pool_order[0].second.cpu_second * 100 << "%" << std::endl;
diff --git a/source/module_base/vector3.h b/source/module_base/vector3.h
index 5079d7cdbd6..f230f616f3c 100644
--- a/source/module_base/vector3.h
+++ b/source/module_base/vector3.h
@@ -32,6 +32,9 @@ class Vector3
Vector3& operator/=(const T &s) { x/=s; y/=s; z/=s; return *this; }
Vector3 operator -() const { return Vector3(-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& normalize(void){ const T m=norm(); x/=m; y/=m; z/=m; return *this; } // Peize Lin update return 2019-09-08
diff --git a/source/module_cell/atom_spec.cpp b/source/module_cell/atom_spec.cpp
index b44960309a0..e0c7565e910 100644
--- a/source/module_cell/atom_spec.cpp
+++ b/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[na];
taud = new ModuleBase::Vector3[na];
- vel = new ModuleBase::Vector3[na];
+ vel = new ModuleBase::Vector3[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]);
diff --git a/source/module_cell/atom_spec.h b/source/module_cell/atom_spec.h
index 75148bb540c..a0935c4c0ed 100644
--- a/source/module_cell/atom_spec.h
+++ b/source/module_cell/atom_spec.h
@@ -41,7 +41,7 @@ class Atom: public Atom_pseudo
ModuleBase::Vector3 *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);
diff --git a/source/module_cell/pseudo_nc.cpp b/source/module_cell/pseudo_nc.cpp
index 162c6ee55e1..bbc5e8270e2 100644
--- a/source/module_cell/pseudo_nc.cpp
+++ b/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 ----------------------";
}
diff --git a/source/module_cell/pseudo_nc.h b/source/module_cell/pseudo_nc.h
index 59af9aeeeac..29fa0d81af0 100644
--- a/source/module_cell/pseudo_nc.h
+++ b/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);
};
diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp
index b1f1b85ac89..ba111be27c9 100644
--- a/source/module_cell/read_atoms.cpp
+++ b/source/module_cell/read_atoms.cpp
@@ -323,7 +323,7 @@ int UnitCell_pseudo::read_atom_species(std::ifstream &ifa, std::ofstream &ofs_ru
latvec.e22 = bba * sinab;
latvec.e23 = 0.0;
latvec.e31 = cba * cosac;
- latvec.e32 = cba * (cosbc - cosac*cosab/sinab);
+ latvec.e32 = cba * (cosbc - cosac*cosab) / sinab;
term = 1.0 + 2.0 * cosab*cosac*cosbc - cosab*cosab - cosac*cosac - cosbc*cosbc;
term = sqrt(term)/sinab;
latvec.e33 = cba * term;
@@ -618,7 +618,7 @@ bool UnitCell_pseudo::read_atom_positions(std::ifstream &ifpos, std::ofstream &o
mv.z = true ;
atoms[it].vel[ia].set(0,0,0);
#ifndef __CMD
- atoms[it].mag[ia]=magnet.start_magnetization[it];
+ //atoms[it].mag[ia]=magnet.start_magnetization[it];//if this line is used, default startmag_type would be 2
#endif
atoms[it].angle1[ia]=0;
atoms[it].angle2[ia]=0;
@@ -925,9 +925,9 @@ bool UnitCell_pseudo::check_tau(void)const
}
#ifdef __LCAO
-void UnitCell_pseudo::print_stru_file(const LCAO_Orbitals &orb, const std::string &fn, const int &type)const
+void UnitCell_pseudo::print_stru_file(const LCAO_Orbitals &orb, const std::string &fn, const int &type, const int &level)const
#else
-void UnitCell_pseudo::print_stru_file(const std::string &fn, const int &type)const
+void UnitCell_pseudo::print_stru_file(const std::string &fn, const int &type, const int &level)const
#endif
{
ModuleBase::TITLE("UnitCell_pseudo","print_stru_file");
@@ -987,19 +987,29 @@ void UnitCell_pseudo::print_stru_file(const std::string &fn, const int &type)con
ofs << "0" << " #magnetism" << std::endl;
#endif
- //2015-05-07, modify
- //ofs << atoms[it].nwl << " #max angular momentum" << std::endl;
- //xiaohui modify 2015-03-15
- //for(int l=0; l<=atoms[it].nwl; l++)
- //{
- // ofs << atoms[it].l_nchi[l] << " #number of zeta for l=" << l << std::endl;
- //}
ofs << atoms[it].na << " #number of atoms" << std::endl;
for(int ia=0; iaprint_cell(ofs, outp);
+ this->print_cell(ofs);
for (int i = 0;i < ntype;i++)
{
- atoms[i].print_Atom(ofs, outp);
+ atoms[i].print_Atom(ofs);
}
ofs.close();
diff --git a/source/module_cell/setup_nonlocal.cpp b/source/module_cell/setup_nonlocal.cpp
index e0432dc318e..58f6926a2e4 100644
--- a/source/module_cell/setup_nonlocal.cpp
+++ b/source/module_cell/setup_nonlocal.cpp
@@ -79,7 +79,7 @@ void InfoNonlocal::Set_NonLocal(
dk,
dr_uniform); // delta k mesh in reciprocal space
- tmpBeta_lm[p1].plot(GlobalV::MY_RANK);
+ if(GlobalV::out_element_info)tmpBeta_lm[p1].plot(GlobalV::MY_RANK);
delete[] beta_r;
@@ -189,7 +189,7 @@ void InfoNonlocal::Set_NonLocal(
dk,
dr_uniform); // delta k mesh in reciprocal space
- tmpBeta_lm[p1].plot(GlobalV::MY_RANK);
+ if(GlobalV::out_element_info)tmpBeta_lm[p1].plot(GlobalV::MY_RANK);
delete[] beta_r;
}
@@ -443,7 +443,7 @@ void InfoNonlocal::Read_NonLocal(
dk,
dr_uniform); // delta k mesh in reciprocal space
- tmpBeta_lm[p1].plot(my_rank);
+ if(GlobalV::out_element_info)tmpBeta_lm[p1].plot(my_rank);
delete[] radial_ps;
delete[] rab_ps;
diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp
index 21a1a2a6e14..9e5170137d0 100644
--- a/source/module_cell/unitcell.cpp
+++ b/source/module_cell/unitcell.cpp
@@ -148,7 +148,7 @@ void UnitCell::bcast_unitcell2(void)
}
#endif
-void UnitCell::print_cell(std::ofstream &ofs, output &outp)const
+void UnitCell::print_cell(std::ofstream &ofs)const
{
if (GlobalV::test_unitcell) ModuleBase::TITLE("UnitCell","print_cell");
@@ -162,10 +162,10 @@ void UnitCell::print_cell(std::ofstream &ofs, output &outp)const
ModuleBase::GlobalFunc::OUT(ofs,"tpiba",tpiba);
ModuleBase::GlobalFunc::OUT(ofs,"omega",omega);
- outp.printM3(ofs,"Lattices Vector (R) : ", latvec);
- outp.printM3(ofs ,"Supercell lattice vector : ", latvec_supercell);
- outp.printM3(ofs, "Reciprocal lattice Vector (G): ", G);
- outp.printM3(ofs, "GGT : ", GGT);
+ output::printM3(ofs,"Lattices Vector (R) : ", latvec);
+ output::printM3(ofs ,"Supercell lattice vector : ", latvec_supercell);
+ output::printM3(ofs, "Reciprocal lattice Vector (G): ", G);
+ output::printM3(ofs, "GGT : ", GGT);
ofs<* posd_in)
+{
+ int iat = 0;
+ for(int it = 0; it < this->ntype; ++it)
+ {
+ Atom* atom = &this->atoms[it];
+ for(int ia = 0; ia < atom->na; ++ia)
+ {
+ if(atom->mbl[ia].x!=0)
+ {
+ atom->tau[ia].x = posd_in[iat].x / this->lat0;
+ }
+ if(atom->mbl[ia].y!=0)
+ {
+ atom->tau[ia].y = posd_in[iat].y / this->lat0;
+ }
+ if(atom->mbl[ia].z!=0)
+ {
+ atom->tau[ia].z = posd_in[iat].z / this->lat0;
+ }
+
+ // the direct coordinates also need to be updated.
+ atom->taud[ia] = atom->tau[ia] * this->GT;
+ iat++;
+ }
+ }
+ assert(iat == this->nat);
+ return;
+}
+
void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in)
{
int iat = 0;
@@ -318,6 +348,21 @@ void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in)
this->periodic_boundary_adjustment();
}
+void UnitCell::update_vel(const ModuleBase::Vector3* vel_in)
+{
+ int iat = 0;
+ for(int it=0; itntype; ++it)
+ {
+ Atom* atom = &this->atoms[it];
+ for(int ia=0; iana; ++ia)
+ {
+ this->atoms[it].vel[ia] = vel_in[iat];
+ ++iat;
+ }
+ }
+ assert(iat == this->nat);
+}
+
void UnitCell::periodic_boundary_adjustment()
{
//----------------------------------------------
@@ -386,7 +431,7 @@ void UnitCell::save_cartesian_position(double* pos)const
return;
}
-bool UnitCell::judge_big_cell(void)
+bool UnitCell::judge_big_cell(void)const
{
double diameter = 2*GlobalV::SEARCH_RADIUS;
diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h
index 5fbc65b0a85..e0532211ffc 100644
--- a/source/module_cell/unitcell.h
+++ b/source/module_cell/unitcell.h
@@ -74,17 +74,19 @@ class UnitCell
public:
UnitCell();
~UnitCell();
- void print_cell(std::ofstream &ofs, output &outp)const;
+ void print_cell(std::ofstream &ofs)const;
void print_cell_xyz(const std::string &fn)const;
void print_cell_cif(const std::string &fn)const;
void update_pos_tau(const double* pos);
+ void update_pos_tau(const ModuleBase::Vector3* posd_in);
void update_pos_taud(const ModuleBase::Vector3* posd_in);
+ void update_vel(const ModuleBase::Vector3* vel_in);
void periodic_boundary_adjustment();
void bcast_atoms_tau();
void save_cartesian_position(double* pos)const;
- bool judge_big_cell(void);
+ bool judge_big_cell(void)const;
protected:
diff --git a/source/module_cell/unitcell_pseudo.cpp b/source/module_cell/unitcell_pseudo.cpp
index 70db502979c..640b327b7bb 100644
--- a/source/module_cell/unitcell_pseudo.cpp
+++ b/source/module_cell/unitcell_pseudo.cpp
@@ -28,7 +28,6 @@ void UnitCell_pseudo::setup_cell(
LCAO_Orbitals &orb,
#endif
const std::string &s_pseudopot_dir,
- output &outp,
const std::string &fn,
std::ofstream &log)
{
@@ -99,7 +98,7 @@ void UnitCell_pseudo::setup_cell(
ok2 = this->read_atom_positions(ifa, log, GlobalV::ofs_warning);
#endif
- if(ok2)
+ if(ok2&&GlobalV::out_element_info)
{
for(int i=0;intype;i++)
{
@@ -169,8 +168,8 @@ void UnitCell_pseudo::setup_cell(
this->invGGT0 = GGT.Inverse();
log << std::endl;
- outp.printM3(log,"Lattice vectors: (Cartesian coordinate: in unit of a_0)",latvec);
- outp.printM3(log,"Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)",G);
+ output::printM3(log,"Lattice vectors: (Cartesian coordinate: in unit of a_0)",latvec);
+ output::printM3(log,"Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)",G);
// OUT(log,"lattice center x",latcenter.x);
// OUT(log,"lattice center y",latcenter.y);
// OUT(log,"lattice center z",latcenter.z);
@@ -197,7 +196,7 @@ void UnitCell_pseudo::setup_cell(
this->read_cell_pseudopots(s_pseudopot_dir, log);
- if(GlobalV::MY_RANK == 0)
+ if(GlobalV::MY_RANK == 0 && GlobalV::out_element_info)
{
for(int it=0; itntype; it++)
{
@@ -380,7 +379,7 @@ void UnitCell_pseudo::setup_cell_classic(
#else
ok2 = this->read_atom_positions(ifa, ofs_running, ofs_warning);
#endif
- if(ok2)
+ if(ok2&&GlobalV::out_element_info)
{
for(int i=0;intype;i++)
{
@@ -422,6 +421,21 @@ void UnitCell_pseudo::setup_cell_classic(
ModuleBase::GlobalFunc::OUT(ofs_running,"Volume (A^3)", this->omega * pow(ModuleBase::BOHR_TO_A, 3));
}
+ //==========================================================
+ // Calculate recip. lattice vectors and dot products
+ // latvec have the unit of lat0, but G has the unit 2Pi/lat0
+ //==========================================================
+ this->GT = latvec.Inverse();
+ this->G = GT.Transpose();
+ this->GGT = G * GT;
+ this->invGGT = GGT.Inverse();
+
+ //LiuXh add 20180515
+ this->GT0 = latvec.Inverse();
+ this->G0 = GT.Transpose();
+ this->GGT0 = G * GT;
+ this->invGGT0 = GGT.Inverse();
+
this->set_iat2itia();
}
@@ -662,17 +676,8 @@ void UnitCell_pseudo::cal_natomwfc(std::ofstream &log)
//LiuXh add a new function here,
//20180515
-void UnitCell_pseudo::setup_cell_after_vc(
- const std::string &s_pseudopot_dir,
- output &outp,
- const std::string &fn, std::ofstream &log)
+void UnitCell_pseudo::setup_cell_after_vc(std::ofstream &log)
{
- if(GlobalV::MY_RANK == 0)
- {
- //std::ifstream ifa(fn.c_str(), ios::in);
- //this->read_atom_species_after_vc(ifa);
- }
-
assert(lat0 > 0.0);
this->omega = abs(latvec.Det()) * this->lat0 * lat0 * lat0;
if(this->omega <= 0)
@@ -722,8 +727,8 @@ Parallel_Common::bcast_double( atom->taud[ia].z );
#endif
log << std::endl;
- outp.printM3(log,"Lattice vectors: (Cartesian coordinate: in unit of a_0)",latvec);
- outp.printM3(log,"Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)",G);
+ output::printM3(log,"Lattice vectors: (Cartesian coordinate: in unit of a_0)",latvec);
+ output::printM3(log,"Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)",G);
return;
}
diff --git a/source/module_cell/unitcell_pseudo.h b/source/module_cell/unitcell_pseudo.h
index d9f1353ff7a..99850dc379f 100644
--- a/source/module_cell/unitcell_pseudo.h
+++ b/source/module_cell/unitcell_pseudo.h
@@ -40,7 +40,6 @@ class UnitCell_pseudo : public UnitCell
LCAO_Orbitals &orb,
#endif
const std::string &s_pseudopot_dir,
- output &outp,
const std::string &fn,
std::ofstream &log);
void setup_cell_classic(
@@ -62,12 +61,12 @@ class UnitCell_pseudo : public UnitCell
int find_type(const std::string &label);
void print_tau(void)const;
#ifdef __LCAO
- void print_stru_file(const LCAO_Orbitals &orb, const std::string &fn, const int &type=1)const; // mohan add 2011-03-22
+ void print_stru_file(const LCAO_Orbitals &orb, const std::string &fn, const int &type=1, const int &level=0)const; // mohan add 2011-03-22
#else
- void print_stru_file(const std::string &fn, const int &type=1)const; // mohan add 2011-03-22
+ void print_stru_file(const std::string &fn, const int &type=1, const int &level=0)const; // mohan add 2011-03-22
#endif
void check_dtau(void);
- void setup_cell_after_vc(const std::string &s_pseudopot_dir, output &outp, const std::string &fn, std::ofstream &log); //LiuXh add 20180515
+ void setup_cell_after_vc(std::ofstream &log); //LiuXh add 20180515
bool set_atom_flag;//added on 2009-3-8 by mohan
@@ -84,7 +83,7 @@ class UnitCell_pseudo : public UnitCell
//void cal_nelec();
void cal_meshx();
void cal_natomwfc(std::ofstream &log);
- void print_unitcell_pseudo(const std::string &fn, output &outp);
+ void print_unitcell_pseudo(const std::string &fn);
bool check_tau(void)const; //mohan add 2011-03-03
bool if_atoms_can_move()const;
bool if_cell_can_change()const;
diff --git a/source/module_md/CMakeLists.txt b/source/module_md/CMakeLists.txt
index 473cb229a92..39f70224d47 100644
--- a/source/module_md/CMakeLists.txt
+++ b/source/module_md/CMakeLists.txt
@@ -1,12 +1,16 @@
add_library(
md
OBJECT
- MD_basic.cpp
- MD_fire.cpp
MD_func.cpp
- MD_thermo.cpp
run_md_classic.cpp
LJ_potential.cpp
DP_potential.cpp
cmd_neighbor.cpp
+ verlet.cpp
+ NVE.cpp
+ MSST.cpp
+ NVT_ADS.cpp
+ NVT_NHC.cpp
+ FIRE.cpp
+ Langevin.cpp
)
diff --git a/source/module_md/DP_potential.cpp b/source/module_md/DP_potential.cpp
index adcb8e36c49..95d77fc7db9 100644
--- a/source/module_md/DP_potential.cpp
+++ b/source/module_md/DP_potential.cpp
@@ -8,7 +8,7 @@ DP_potential::DP_potential(){}
DP_potential::~DP_potential(){}
-void DP_potential::DP_pot(UnitCell_pseudo &ucell_c, double &potential, ModuleBase::Vector3 *force, ModuleBase::matrix &stress)
+void DP_potential::DP_pot(const UnitCell_pseudo &ucell_c, double &potential, ModuleBase::Vector3 *force, ModuleBase::matrix &stress)
{
#ifdef __DPMD
ModuleBase::TITLE("DP_potential", "DP_pot");
@@ -57,7 +57,7 @@ void DP_potential::DP_pot(UnitCell_pseudo &ucell_c, double &potential, ModuleBas
#endif
}
-void DP_potential::DP_init(UnitCell_pseudo &ucell_c,
+void DP_potential::DP_init(const UnitCell_pseudo &ucell_c,
std::vector &cell,
std::vector &coord,
std::vector &atype)
diff --git a/source/module_md/DP_potential.h b/source/module_md/DP_potential.h
index 9fa5aa3e291..ec2b605df81 100644
--- a/source/module_md/DP_potential.h
+++ b/source/module_md/DP_potential.h
@@ -12,12 +12,12 @@ class DP_potential
DP_potential();
~DP_potential();
- static void DP_init(UnitCell_pseudo &ucell_c,
+ static void DP_init(const UnitCell_pseudo &ucell_c,
std::vector &cell,
std::vector &coord,
std::vector &atype);
- static void DP_pot(UnitCell_pseudo &ucell_c,
+ static void DP_pot(const UnitCell_pseudo &ucell_c,
double &potential,
ModuleBase::Vector3 *force,
ModuleBase::matrix &stress);
diff --git a/source/module_md/FIRE.cpp b/source/module_md/FIRE.cpp
new file mode 100644
index 00000000000..8a9db220076
--- /dev/null
+++ b/source/module_md/FIRE.cpp
@@ -0,0 +1,209 @@
+#include "FIRE.h"
+#include "MD_func.h"
+
+FIRE::FIRE(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para_in, unit_in)
+{
+ dt_max = -1.0;
+ alpha_start = 0.10;
+ alpha = alpha_start;
+
+ finc = 1.1;
+ fdec = 0.5;
+ f_alpha = 0.99;
+ N_min = 4;
+ negative_count = 0;
+}
+
+FIRE::~FIRE(){}
+
+void FIRE::setup()
+{
+ ModuleBase::TITLE("FIRE", "setup");
+ ModuleBase::timer::tick("FIRE", "setup");
+
+ Verlet::setup();
+
+ check_force();
+
+ ModuleBase::timer::tick("FIRE", "setup");
+}
+
+void FIRE::first_half()
+{
+ ModuleBase::TITLE("FIRE", "first_half");
+ ModuleBase::timer::tick("FIRE", "first_half");
+
+ for(int i=0; i> step_rst_ >> alpha >> negative_count >> dt_max >> mdp.dt;
+
+ file.close();
+ }
+
+#ifdef __MPI
+ MPI_Bcast(&step_rst_, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&alpha, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&negative_count, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&dt_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&mdp.dt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+#endif
+}
+
+void FIRE::check_force()
+{
+ max = 0;
+
+ for(int i=0; i 0)
+ {
+ negative_count++;
+ if(negative_count >= N_min)
+ {
+ mdp.dt = min(mdp.dt*finc, dt_max);
+ alpha *= f_alpha;
+ }
+ }
+ else
+ {
+ mdp.dt *= fdec;
+ negative_count = 0;
+
+ for(int i=0; i *force,
ModuleBase::matrix &stress)
@@ -34,7 +34,7 @@ double LJ_potential::Lennard_Jones(UnitCell_pseudo &ucell_c,
distance = dtau.norm();
if(distance <= INPUT.mdp.rcut_lj)
{
- potential += LJ_energy(distance) - LJ_energy(INPUT.mdp.rcut_lj);
+ potential += LJ_energy(distance); // - LJ_energy(INPUT.mdp.rcut_lj);
ModuleBase::Vector3 f_ij = LJ_force(distance, dtau);
force[index] += f_ij;
LJ_virial(virial, f_ij, dtau);
@@ -56,7 +56,7 @@ double LJ_potential::Lennard_Jones(UnitCell_pseudo &ucell_c,
return potential/2.0;
}
-double LJ_potential::Lennard_Jones(UnitCell_pseudo &ucell_c,
+double LJ_potential::Lennard_Jones(const UnitCell_pseudo &ucell_c,
CMD_neighbor &cmd_neigh,
ModuleBase::Vector3 *force,
ModuleBase::matrix &stress)
@@ -95,7 +95,7 @@ double LJ_potential::Lennard_Jones(UnitCell_pseudo &ucell_c,
if(distance <= INPUT.mdp.rcut_lj)
{
- potential += LJ_energy(distance) - LJ_energy(INPUT.mdp.rcut_lj);
+ potential += LJ_energy(distance); // - LJ_energy(INPUT.mdp.rcut_lj);
ModuleBase::Vector3 f_ij = LJ_force(distance, temp*ucell_c.lat0);
force[i] = force[i] + f_ij;
LJ_virial(virial, f_ij, temp*ucell_c.lat0);
diff --git a/source/module_md/LJ_potential.h b/source/module_md/LJ_potential.h
index 3cc8c4da1bd..8a0055caea9 100644
--- a/source/module_md/LJ_potential.h
+++ b/source/module_md/LJ_potential.h
@@ -13,8 +13,8 @@ class LJ_potential
LJ_potential();
~LJ_potential();
- static double Lennard_Jones(UnitCell_pseudo &ucell_c, Grid_Driver &grid_neigh, ModuleBase::Vector3 *force, ModuleBase::matrix &stress);
- static double Lennard_Jones(UnitCell_pseudo &ucell_c, CMD_neighbor &cmd_neigh, ModuleBase::Vector3 *force, ModuleBase::matrix &stress);
+ static double Lennard_Jones(const UnitCell_pseudo &ucell_c, Grid_Driver &grid_neigh, ModuleBase::Vector3 *force, ModuleBase::matrix &stress);
+ static double Lennard_Jones(const UnitCell_pseudo &ucell_c, CMD_neighbor &cmd_neigh, ModuleBase::Vector3 *force, ModuleBase::matrix &stress);
static double LJ_energy(const double d);
static ModuleBase::Vector3 LJ_force(const double d, const ModuleBase::Vector3 dr);
static void LJ_virial(double *virial, const ModuleBase::Vector3 &force, const ModuleBase::Vector3 &dtau);
diff --git a/source/module_md/Langevin.cpp b/source/module_md/Langevin.cpp
new file mode 100644
index 00000000000..2ff884ab364
--- /dev/null
+++ b/source/module_md/Langevin.cpp
@@ -0,0 +1,79 @@
+#include "Langevin.h"
+#include "MD_func.h"
+
+Langevin::Langevin(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para_in, unit_in)
+{
+ // convert to a.u. unit
+ mdp.damp /= ModuleBase::AU_to_FS;
+}
+
+Langevin::~Langevin(){}
+
+void Langevin::setup()
+{
+ ModuleBase::TITLE("Langevin", "setup");
+ ModuleBase::timer::tick("Langevin", "setup");
+
+ Verlet::setup();
+
+ post_force();
+
+ ModuleBase::timer::tick("Langevin", "setup");
+}
+
+void Langevin::first_half()
+{
+ ModuleBase::TITLE("Langevin", "first_half");
+ ModuleBase::timer::tick("Langevin", "first_half");
+
+ Verlet::first_half();
+
+ ModuleBase::timer::tick("Langevin", "first_half");
+}
+
+void Langevin::second_half()
+{
+ ModuleBase::TITLE("Langevin", "second_half");
+ ModuleBase::timer::tick("Langevin", "second_half");
+
+ post_force();
+
+ Verlet::second_half();
+
+ ModuleBase::timer::tick("Langevin", "second_half");
+}
+
+void Langevin::outputMD()
+{
+ Verlet::outputMD();
+}
+
+void Langevin::write_restart()
+{
+ Verlet::write_restart();
+}
+
+void Langevin::restart()
+{
+ Verlet::restart();
+}
+
+void Langevin::post_force()
+{
+ temp_target();
+
+ for(int i=0; i[ucell.nat];
- cart_change=new ModuleBase::Vector3[ucell.nat];
- tauDirectChange=new ModuleBase::Vector3[ucell.nat];
- allmass=new double[ucell.nat];
- ionmbl=new ModuleBase::Vector3[ucell.nat];
-
- frozen_freedom_ = mdf.getMassMbl(ucell, mdp, allmass, ionmbl);
- mdf.InitVel(ucell, temperature_, allmass, frozen_freedom_, ionmbl, vel);
- // if (ucell.set_vel) // Yuanbo Li 2021-08-01
- // {
- // int iat=0; //initialize velocity of atoms from STRU liuyu 2021-07-14
- // for(int it=0; it *force, const ModuleBase::matrix &stress)
-{
-//------------------------------------------------------------------------------
-// DESCRIPTION:
-// Molecular dynamics calculation with fixed Volume and slight fluctuated temperature
-// Using thermostat : 1, Nose-Hoover Chains; 2, Langevin; 3, Anderson
-// Normal Nose-Hoover thermostat method is retained for test.
-//------------------------------------------------------------------------------
-
- ModuleBase::TITLE("MD_basic","runnvt");
- ModuleBase::timer::tick("MD_basic","runnvt");
- step_=step1+step_rst_;
- //the real MD step
-
- //change target temperature
-// if (step_!=1)mdf.ReadNewTemp( step_ );
-
- std::cout << " ------------------------------------------------------------" << std::endl;
- std::cout << " Target temperature : " << temperature_*ModuleBase::Hartree_to_K << " (K)"<< std::endl;
-
- if(step_==1||step_%mdp.fixTemperature==1)
- {
- for(int k=0;kupdate_half_velocity(force);
- }
- else
- {
- this->update_half_velocity(force);
-
- twiceKE=mdf.GetAtomKE(ucell.nat, vel, allmass);
- twiceKE = 2 * twiceKE;
- if(mdp.NVT_control==1)
- {
- hamiltonian = mdt.NHChamiltonian(
- twiceKE/2,
- energy_,
- temperature_,
- frozen_freedom_);
- }
- else hamiltonian = mdf.Conserved(twiceKE/2, energy_, ucell.nat-frozen_freedom_);
- //Note: side scheme position before
- //Now turn to middle scheme.
- this->update_half_velocity(force);
- }
-
- // Update the Non-Wrapped cartesion coordinates
- if(mdp.mdtype==2) mdf.scalevel(ucell.nat, frozen_freedom_, temperature_, vel, allmass);//choose velocity scaling method
-
- update_half_direct(1);
-
- mdt.Integrator(mdp.NVT_control, temperature_, vel, allmass);//thermostat interact with velocity
- twiceKE=mdf.GetAtomKE(ucell.nat, vel, allmass);
- twiceKE = 2 * twiceKE;
-
- update_half_direct(0);
-
- // Calculate the maximal velocities.
- // The step in fractional coordinates
- double maxStep = 0;
- for(int ii=0;iimaxStep)
- {
- maxStep = pow(vel[ii].x,2)+pow(vel[ii].y,2)+pow(vel[ii].z,2);
- }
- }
- getTaudUpdate();
-
- //save the atom position change to DFT module
- save_output_position();
- maxStep = sqrt(maxStep)*mdp.dt;
-
- if (!GlobalV::MY_RANK)
- {
- GlobalV::ofs_running << " maxForce : " << std::setw(10) << maxForce << std::endl;
- GlobalV::ofs_running << " maxStep : " << std::setw(10) << maxStep << std::endl;
- GlobalV::ofs_running << " " < *force, const ModuleBase::matrix &stress)
-{
-//-------------------------------------------------------------------------------
-// Adiabatic ensemble
-// Molecular dynamics calculation with Verlet algorithm
-//-------------------------------------------------------------------------------
-
- ModuleBase::TITLE("MD_basic","runNVE");
- ModuleBase::timer::tick("MD_basic","runNVE");
- step_=step1+step_rst_;
-
- // Set up extra output to ion optimizer / MD header
- //std::cout<<"Time interval : "<< mdp.dt*fundamentalTime*1e15<< " (fs)"<update_half_velocity(force);
- // (2) 2nd step of Verlet-Velocity
- // Update the Non-Wrapped cartesion coordinate
- twiceKE = mdf.GetAtomKE(ucell.nat, vel, allmass);
- twiceKE = 2 * twiceKE;
- if(step_!=1||mdp.rstMD==1)this->update_half_velocity(force);
- update_half_direct(1);
- update_half_direct(0);
- // Calculate the maximal velocities.
- // The step in fractional coordinates
- double maxStep = 0;
- for(int i = 0;i< ucell.nat;i++)
- {
- if((pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0))>maxStep)
- {
- maxStep = pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0);
- }
- }
- getTaudUpdate();
- save_output_position();
- maxStep = sqrt(maxStep)*mdp.dt;
-
-
- // calculate the conserved quantity during MD
- double hamiltonian = mdf.Conserved(twiceKE/2, energy_, 3*ucell.nat-frozen_freedom_);
-
- //std::cout<< std::setprecision (9)< *force, const ModuleBase::matrix &stress)
-{
-//-------------------------------------------------------------------------------
-// REFERENCES:
-//
- ModuleBase::TITLE("MD_basic","runFIRE");
- ModuleBase::timer::tick("MD_basic","runFIRE");
- step_ = step1;
- // if(step_==1)
- // {
- // std::cout<<" FIRE Start."<update_half_velocity(force);
-
- // (2) 2nd step of Verlet-Velocity
- // Update the Non-Wrapped cartesion coordinate
- twiceKE = mdf.GetAtomKE(ucell.nat, vel, allmass);
- twiceKE = 2 * twiceKE;
- if(step_!=1)this->update_half_velocity(force);
-
- double largest_grad_FIRE = 0.0;
- for(int i=0;imaxStep)
- {
- maxStep = pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0);
- }
- }
- getTaudUpdate();
-
- save_output_position();
- maxStep = sqrt(maxStep)*mdp.dt;
-
- double hamiltonian = mdf.Conserved(twiceKE/2, energy_, 3 * ucell.nat - frozen_freedom_);
-
- // Output the message to the screen.
- if (!GlobalV::MY_RANK)
- {
- GlobalV::ofs_running << " ------------------------------------------------------------" << std::endl;
- GlobalV::ofs_running << " "< *force)
-{
- for(int ii=0;iistep_;
-}
-
-//output pressure of total MD system, P = tr(stress) + P_kin
-void MD_basic::outStressMD(const ModuleBase::matrix& stress, const double& twiceKE)
-{
- GlobalV::ofs_running<<"\noutput Pressure for check!"< fracStep;
- for(int ii=0;ii *force, const ModuleBase::matrix &stress);//NVT ensemble MD
- void runNVE(int step1, double potential, ModuleBase::Vector3 *force, const ModuleBase::matrix &stress); //NVE ensemble MD
- bool runFIRE(int step1, double potential, ModuleBase::Vector3 *force, const ModuleBase::matrix &stress); //relax method FIRE
- int getRealStep();
-
- private:
- MD_func mdf;
- MD_thermo mdt;
- MD_parameters &mdp;
- UnitCell_pseudo &ucell;
- MD_fire fire;
-
- double temperature_;
- const static double fundamentalTime;
- int outputstressperiod_;
- int step_rst_;
- int step_;
- double energy_;
- int frozen_freedom_;
- double oldEtot_;
-
- //allocable variaints
- ModuleBase::Vector3 *vel;// velocity of each atom, unit is a.u.
- double *allmass; //mass of each atom
- //ModuleBase::Vector3 *force; //force of each atom
- //matrix stress; //stress for this lattice
- ModuleBase::Vector3 *ionmbl; //if frozen of each atom
- ModuleBase::Vector3 *tauDirectChange; //change of dirac coord of atoms
- ModuleBase::Vector3 *cart_change; //cartensian coord of atoms, *not* wrapped
-
- //repete code
- void update_half_velocity(ModuleBase::Vector3 *force);
- void update_half_direct(const bool is_restart);
- void save_output_position();
- void outStressMD(const ModuleBase::matrix& stress, const double& twiceKE);
- void getTaudUpdate();
-};
-
-#endif
diff --git a/source/module_md/MD_fire.cpp b/source/module_md/MD_fire.cpp
deleted file mode 100644
index 4b20346556a..00000000000
--- a/source/module_md/MD_fire.cpp
+++ /dev/null
@@ -1,69 +0,0 @@
-#include "MD_fire.h"
-
-MD_fire::MD_fire()
-{
- dt_max=-1.0;
- alpha_start=0.10;
- alpha = alpha_start;
-
- finc=1.1;
- fdec=0.5;
- f_alpha=0.99;
- N_min=4;
- negative_count=0;
-}
-
-void MD_fire::check_FIRE(
- const int& numIon,
- const ModuleBase::Vector3* force,
- double& deltaT,
- ModuleBase::Vector3* vel)
-{
- if(dt_max<0) dt_max = deltaT * 2.5; //initial dt_max
- double P(0.0);
- double sumforce(0.0);
- double normvel(0.0);
-
- for(int iatom =0;iatom 0 )
- {
- negative_count +=1 ;
- if(negative_count >=N_min)
- {
- deltaT=min(deltaT*finc,dt_max);
- alpha= alpha *f_alpha;
- }
- }
- else if( P<=0)
- {
- deltaT=deltaT*fdec;
- negative_count = 0;
-
- for(int iatom=0; iatom* force,
- double& deltaT,
- ModuleBase::Vector3* vel);
- private:
- double alpha_start ;//alpha_start begin
- double alpha;//alpha begin
- double finc;//finc begin
- double fdec;//fdec begin
- double f_alpha;
- int N_min;
- double dt_max;//dt_max
- int negative_count;//Negative count
-};
-
-#endif
diff --git a/source/module_md/MD_func.cpp b/source/module_md/MD_func.cpp
index 6c6ec2b43f1..61e26ab7908 100644
--- a/source/module_md/MD_func.cpp
+++ b/source/module_md/MD_func.cpp
@@ -1,4 +1,50 @@
#include "MD_func.h"
+#include "cmd_neighbor.h"
+#include "LJ_potential.h"
+#include "DP_potential.h"
+#include "../input.h"
+#include "../module_neighbor/sltk_atom_arrange.h"
+#include "../module_neighbor/sltk_grid_driver.h"
+#include "../module_base/global_variable.h"
+
+#ifndef __CMD
+#include "../src_pw/run_md_pw.h"
+#endif
+
+#ifdef __LCAO
+#include "../src_lcao/run_md_lcao.h"
+#endif
+
+
+double MD_func::gaussrand()
+{
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+
+ if ( phase == 0 )
+ {
+ do
+ {
+ double U1 = rand()/double(RAND_MAX);
+ double U2 = rand()/double(RAND_MAX);
+
+ V1 = 2.0 * U1 - 1.0;
+ V2 = 2.0 * U2 - 1.0;
+ S = V1 * V1 + V2 * V2;
+ } while(S >= 1 || S == 0);
+
+ X = V1 * sqrt(-2.0 * log(S) / S);
+ }
+ else
+ {
+ X = V2 * sqrt(-2.0 * log(S) / S);
+ }
+
+ phase = 1 - phase;
+
+ return X;
+}
bool MD_func::RestartMD(const int& numIon, ModuleBase::Vector3* vel, int& step_rst)
{
@@ -94,20 +140,64 @@ void MD_func::mdRestartOut(const int& step, const int& recordFreq, const int& nu
return;
}
-double MD_func::GetAtomKE(const int& numIon, const ModuleBase::Vector3* vel, const double * allmass){
+double MD_func::GetAtomKE(
+ const int &numIon,
+ const ModuleBase::Vector3 *vel,
+ const double *allmass)
+{
+ double ke = 0;
+
+ for(int ion=0; ion *vel,
+ const double *allmass,
+ double &kinetic,
+ ModuleBase::matrix &stress)
+{
//---------------------------------------------------------------------------
// DESCRIPTION:
-// This function calculates the classical kinetic energy of a group of atoms.
+// This function calculates the classical kinetic energy of atoms
+// and its contribution to stress.
//----------------------------------------------------------------------------
- double ke = 0.0;
+ kinetic = MD_func::GetAtomKE(unit_in.nat, vel, allmass);
+
+ ModuleBase::matrix temp;
+ temp.create(3,3); // initialize
- // kinetic energy = \sum_{i} 1/2*m_{i}*(vx^2+vy^2+vz^2)
- for(int ion=0;ion *force,
+ ModuleBase::matrix &stress)
+{
+ ModuleBase::TITLE("MD_func", "force_stress");
+ ModuleBase::timer::tick("MD_func", "force_stress");
+
+ if(mdp.md_potential == "LJ")
+ {
+ bool which_method = unit_in.judge_big_cell();
+ if(which_method)
+ {
+ CMD_neighbor cmd_neigh;
+ cmd_neigh.neighbor(unit_in);
+
+ potential = LJ_potential::Lennard_Jones(
+ unit_in,
+ cmd_neigh,
+ force,
+ stress);
}
- for(ion=0;ion *force)
+{
+ if(GlobalV::MY_RANK) return;
+
+ std::stringstream file;
+ file << GlobalV::global_out_dir << "MD_dump";
+ std::ofstream ofs;
+ if(step == 0)
+ {
+ ofs.open(file.str(), ios::trunc);
+ }
+ else
+ {
+ 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;
+
+ ofs << "MDstep: " << step << std::endl;
+
+ ofs << "VIRIAL (KBAR)" << std::endl;
+ ofs << std::setprecision(12) << std::setiosflags(ios::fixed);
+ 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 << "\nFORCE (eV/Angstrom)" << std::endl;
+ for(int i=0; i &frozen,
ModuleBase::Vector3* ionmbl)
{
//some prepared information
//mass and degree of freedom
int ion=0;
- int frozen_freedom=0;
+ frozen.set(0,0,0);
for(int it=0;it* vel, int& step_rst);
- void mdRestartOut(const int& step, const int& recordFreq, const int& numIon, ModuleBase::Vector3* vel);
- double GetAtomKE(const int& numIon, const ModuleBase::Vector3* vel, const double* allmass);
- void InitVelocity(
- const int& numIon,
- const double& temperature,
- const double& fundamentalTime,
- const double* allmass,
- ModuleBase::Vector3* vel);
+ static double gaussrand();
+
+ static bool RestartMD(const int& numIon, ModuleBase::Vector3* vel, int& step_rst);
+ static void mdRestartOut(const int& step, const int& recordFreq, const int& numIon, ModuleBase::Vector3* vel);
- void InitVel(
+ static void InitPos( const UnitCell_pseudo &unit_in, ModuleBase::Vector3* pos);
+
+ static void InitVel(
const UnitCell_pseudo &unit_in,
const double& temperature,
double* allmass,
@@ -29,32 +26,58 @@ class MD_func
ModuleBase::Vector3* ionmbl,
ModuleBase::Vector3* vel);
- void ReadVel(
- const UnitCell_pseudo &unit_in,
- ModuleBase::Vector3* vel);
+ static void ReadVel(const UnitCell_pseudo &unit_in, ModuleBase::Vector3* vel);
- void RandomVel(
+ static void RandomVel(
const int& numIon,
const double& temperature,
const double* allmass,
const int& frozen_freedom,
+ const ModuleBase::Vector3 frozen,
const ModuleBase::Vector3* ionmbl,
ModuleBase::Vector3* vel);
-// void ReadNewTemp(int step);
- std::string intTurnTostring(long int iter,std::string path);
- int getMassMbl(const UnitCell_pseudo &unit_in,
+ static void force_virial(
+ const int &istep,
const MD_parameters &mdp,
+ const UnitCell_pseudo &unit_in,
+ double &potential,
+ ModuleBase::Vector3 *force,
+ ModuleBase::matrix &stress);
+
+ static double GetAtomKE(
+ const int &numIon,
+ const ModuleBase::Vector3 *vel,
+ const double *allmass);
+
+ static void kinetic_stress(
+ const UnitCell_pseudo &unit_in,
+ const ModuleBase::Vector3 *vel,
+ const double *allmass,
+ double &kinetic,
+ ModuleBase::matrix &stress);
+
+ static void outStress(const ModuleBase::matrix &virial, const ModuleBase::matrix &stress);
+
+ static void MDdump(
+ const int &step,
+ const int &natom,
+ const ModuleBase::matrix &virial,
+ const ModuleBase::Vector3 *force);
+
+ static std::string intTurnTostring(long int iter,std::string path);
+ static void getMassMbl(const UnitCell_pseudo &unit_in,
double* allmass,
+ ModuleBase::Vector3 &frozen,
ModuleBase::Vector3* ionmbl);
- void printpos(const std::string& file, const int& iter, const int& recordFreq, const UnitCell_pseudo& unit_in);
- void scalevel(
+ static void printpos(const std::string& file, const int& iter, const int& recordFreq, const UnitCell_pseudo& unit_in);
+ static void scalevel(
const int& numIon,
const int& nfrozen,
const double& temperature,
ModuleBase::Vector3* vel,
const double* allmass);
- double MAXVALF(const int numIon, const ModuleBase::Vector3* force);
- double Conserved(const double KE, const double PE, const int number);
+ static double MAXVALF(const int numIon, const ModuleBase::Vector3* force);
+ static double Conserved(const double KE, const double PE, const int number);
};
#endif
diff --git a/source/module_md/MD_parameters.h b/source/module_md/MD_parameters.h
index a182e0be8b8..337addc2e3e 100644
--- a/source/module_md/MD_parameters.h
+++ b/source/module_md/MD_parameters.h
@@ -8,49 +8,77 @@ class MD_parameters
public:
MD_parameters()
{
- mdtype=1;
+ rstMD = 0;
+ mdtype = 2;
+ dt = 1;
+ tfirst = 0;
+ tlast = -1;
+ recordFreq = 1;
+
+ // useless
NVT_tau=0;
NVT_control=1;
- dt=1;
- Qmass=1;
- tfirst=0; //kelvin
- tlast=-1;
- recordFreq=1;
mdoutputpath="mdoutput";
- rstMD=0;
fixTemperature=1;
ediff=1e-4;
ediffg=1e-3;
- MNHC=4;
+ // Classic MD
md_potential = "FP";
- rcut_lj = 8.5; // \AA
- epsilon_lj = 0.01032; // eV
- sigma_lj = 3.405; // \AA
+ rcut_lj = 8.5;
+ epsilon_lj = 0.01032;
+ sigma_lj = 3.405;
+
+ // MSST
+ direction = 2;
+ Qmass = 1;
+ velocity = 0;
+ viscosity = 0;
+ tscale = 0.01;
+
+ // NHC
+ tfreq = 1;
+ MNHC = 4;
+
+ // Langevin
+ damp = 1;
};
~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
- double Qmass; // Inertia of extended system variable
+ int rstMD; // 1: restart MD, 0: no restart MD
+ int mdtype; // -1: FIRE, 0: NVE, 1: NVT ADS, 2: NVT NHC, 3: NPT, 4: MSST
double dt ; // Time increment (hbar/E_hartree)
- double tfirst; //Temperature (in Hartree, 1 Hartree ~ 3E5 K)
- double tlast;
+ double tfirst; // Temperature (in Hartree, 1 Hartree ~ 3E5 K)
+ double tlast; // Target temperature
int recordFreq; // The period to dump MD information for monitoring and restarting MD
+
+ // useless
int NVT_control;
double NVT_tau;
- int MNHC;
std::string mdoutputpath;// output directory of md files: .ion .vel
double ediff; //parameter for constrain
double ediffg; //parameter for constrain
int fixTemperature;
- /*Classic MD*/ // liuyu 2021-07-30
+ // Classic MD // liuyu 2021-07-30
std::string md_potential; // choose potential: LJ, DP, FP
double rcut_lj; // cutoff radius of LJ potential (\AA)
double epsilon_lj; // the value of epsilon for LJ potential (eV)
double sigma_lj; // the value of sigma for LJ potential (\AA)
+
+ // MSST
+ int direction; // shock direction: 0, 1, 2
+ double velocity; // shock velocity (\AA/fs)
+ double Qmass; // cell mass-like parameter (mass^2/length^4)
+ double viscosity; // artificial viscosity (mass/length/time)
+ double tscale; // reduction in initial temperature (0~1)
+
+ // NHC
+ double tfreq; // Oscillation frequency, used to determine Qmass of NHC
+ int MNHC; // num of NHC
+
+ // Langevin
+ double damp; // damping parameter (time units)
};
diff --git a/source/module_md/MD_thermo.cpp b/source/module_md/MD_thermo.cpp
deleted file mode 100644
index 422ac658ef4..00000000000
--- a/source/module_md/MD_thermo.cpp
+++ /dev/null
@@ -1,609 +0,0 @@
-#include"MD_thermo.h"
-
-MD_thermo::MD_thermo()
-{
- G = new ModuleBase::Vector3[1];
- NHCeta = new ModuleBase::Vector3[1];
- NHCpeta = new ModuleBase::Vector3[1];
-}
-
-MD_thermo::~MD_thermo()
-{
- delete[] G;
- delete[] NHCeta;
- delete[] NHCpeta;
-}
-
-void MD_thermo::init_NHC(
- const int &MNHC_in,
- const double &Qmass_in,
- const double &NVT_tau_in,
- const double &dt_in,
- const int &NVT_control,
- std::ofstream &ofs,
- const int &numIon,
- const double &temperature,
- const ModuleBase::Vector3* vel,
- const double* allmass
- )
-{
- this->MNHC_ = MNHC_in;
- this->Qmass_ = Qmass_in;
- this->NVT_tau_ = NVT_tau_in;
- this->dt_ = dt_in;
- this->numIon_ = numIon;
- unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
- init_by_array(init, length);
- ofs<<" ...............Nose-Hoover Chain parameter initialization............... " << std::endl;
- ofs<<" Temperature = "<< temperature << std::endl;
- ofs<<" Temperature2 = "<< temperature/ModuleBase::K_BOLTZMAN_AU << std::endl;
- ofs<<" NHC frequency = "<< 1.0/NVT_tau_ << std::endl;
- ofs<<" NHC chain = "<< MNHC_ << std::endl;
- ofs<<" Qmass = "<< Qmass_ << std::endl;
- ofs<<" ............................................................... " << std::endl;
- w[0]=0.7845136105;
- w[6]=0.7845136105;
- w[1]=0.2355732134;
- w[5]=0.2355732134;
- w[2]=-1.177679984;
- w[4]=-1.177679984;
- w[3]=1-w[0]-w[1]-w[2]-w[4]-w[5]-w[6];
- //
- for(int i=0;i[MNHC_*numIon_];
- delete[] NHCeta;
- NHCeta=new ModuleBase::Vector3[MNHC_*numIon_];
- delete[] NHCpeta;
- NHCpeta=new ModuleBase::Vector3[MNHC_*numIon_];
-
- for(int j=0;j=1;j--)
- {
- for(int k=0;k= 1 || S == 0);
-
- X = V1 * sqrt(-2.0 * log(S) / S);
- }
- else
- {
- X = V2 * sqrt(-2.0 * log(S) / S);
- }
-
- phase = 1 - phase;
-
- return X;
-}
-
-//zifei
-/* initializes mt[N] with a seed */
-void MD_thermo::init_genrand(unsigned long s)
-{
- mt[0]= s & 0xffffffffUL;
- for (mti=1; mti> 30)) + mti);
- /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
- /* In the previous versions, MSBs of the seed affect */
- /* only MSBs of the array mt[]. */
- /* 2002/01/09 modified by Makoto Matsumoto */
- mt[mti] &= 0xffffffffUL;
- /* for >32 bit machines */
- }
-}
-
-
-/* initialize by an array with array-length */
-/* init_key is the array for initializing keys */
-/* key_length is its length */
-/* slight change for C++, 2004/2/26 */
-void MD_thermo::init_by_array(unsigned long init_key[], int key_length)
-{
- init_genrand(19650218UL);
- int i=1;
- int j=0;
- int k = (N_mt19937>key_length ? N_mt19937 : key_length);
- for (; k; k--)
- {
- mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
- + init_key[j] + j; /* non linear */
- mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
- i++; j++;
- if (i>=N_mt19937) { mt[0] = mt[N_mt19937-1]; i=1; }
- if (j>=key_length) j=0;
- }
- for (k=N_mt19937-1; k; k--)
- {
- mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
- - i; /* non linear */
- mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
- i++;
- if (i>=N_mt19937) { mt[0] = mt[N_mt19937-1]; i=1; }
- }
-
- mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
-}
-
-/* generates a random number on [0,0xffffffff]-interval */
-unsigned long MD_thermo::genrand_int32(void)
-{
- unsigned long y;
- static unsigned long mag01[2]={0x0UL, MATRIX_A};
- /* mag01[x] = x * MATRIX_A for x=0,1 */
-
- if (mti >= N_mt19937)
- { /* generate N words at one time */
- int kk;
-
- if (mti == N_mt19937+1)
- { /* if init_genrand() has not been called, */
- init_genrand(5489UL); /* a default initial seed is used */
- }
- for (kk=0;kk> 1) ^ mag01[y & 0x1UL];
- }
- for (;kk> 1) ^ mag01[y & 0x1UL];
- }
- y = (mt[N_mt19937-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
- mt[N_mt19937-1] = mt[M_mt19937-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
-
- mti = 0;
- }
-
- y = mt[mti++];
-
- /* Tempering */
- y ^= (y >> 11);
- y ^= (y << 7) & 0x9d2c5680UL;
- y ^= (y << 15) & 0xefc60000UL;
- y ^= (y >> 18);
-
- return y;
-}
-
-/* generates a random number on [0,0x7fffffff]-interval */
-long MD_thermo::genrand_int31(void)
-{
- return (long)(genrand_int32()>>1);
-}
-
-/* generates a random number on [0,1]-real-interval */
-double MD_thermo::genrand_real1(void)
-{
- return genrand_int32()*(1.0/4294967295.0);
- /* divided by 2^32-1 */
-}
-
-/* generates a random number on [0,1)-real-interval */
-double MD_thermo::genrand_real2(void)
-{
- return genrand_int32()*(1.0/4294967296.0);
- /* divided by 2^32 */
-}
-
-/* generates a random number on (0,1)-real-interval */
-double MD_thermo::genrand_real3(void)
-{
- return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
- /* divided by 2^32 */
-}
-
-/* generates a random number on [0,1) with 53-bit resolution*/
-double MD_thermo::genrand_res53(void)
-{
- unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
- return(a*67108864.0+b)*(1.0/9007199254740992.0);
-}
-
-
-void MD_thermo::Integrator(
- const int control,
- const double &temperature,
- ModuleBase::Vector3* vel,
- const double* allmass)
-{
- if(control == 1) NHCIntegrator(temperature, vel, allmass);
- else if(control == 2) LGVIntegrator(temperature, vel, allmass);
- else if(control == 3) ADSIntegrator(temperature, vel, allmass);
- else ModuleBase::WARNING_QUIT("MD_thermo:Integrator", "please choose available reservoir!!!");
- return;
-}
-
-void MD_thermo::LGVIntegrator(
- const double &temperature,
- ModuleBase::Vector3* vel,
- const double* allmass
-)
-{
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-// This function propagates the Langevin dynamics
-//--------------------------------------------------------------------
- double randomx,randomy,randomz;
- double tempV;
- double c1k = exp(-1/NVT_tau_ * dt_);
- //c1k=e^(-gamma*dt)
- double c2k=sqrt(1.0-c1k*c1k);
-
- for(int iatom=0;iatom* vel,
- const double* allmass
-)
-{
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-// This function propagates the Andersen thermostat
-//--------------------------------------------------------------------
-
- double uniform_random;
- double ranx;
- double rany;
- double ranz;
-
- double nu = 1/ NVT_tau_;
-
- for(int iatom=0;iatom* vel,
- const double* allmass
-){
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-// This function propagates the Nose-Hoover Chain
-// NHCpeta :NHC momentum
-// NHCeta : NHC position
-// Q : NHC mass
-//--------------------------------------------------------------------
-
- for(int i=0;i=0;j--)
- {
- for(int k=0;k=0;j--)
- {
- for(int k=0;k 1&& step%recordFreq==0 ) )
- pass =1;
- if (!pass) return;
-
- if(!GlobalV::MY_RANK){
- std::stringstream ssc;
- ssc << GlobalV::global_out_dir << "Restart_md.dat";
- std::ofstream file(ssc.str().c_str(), ios::app);
- file<<'\n';
- file<<"MD_THERMOSTAT"<
-void MD_thermo::NHC_restart()
-{
- int error=0;
- double *nhcg=new double[numIon_*3*MNHC_];
- double *nhcp=new double[numIon_*3*MNHC_];
- double *nhce=new double[numIon_*3*MNHC_];
- if (!GlobalV::MY_RANK)
- {
- std::stringstream ssc;
- ssc << GlobalV::global_readin_dir << "Restart_md.dat";
- std::ifstream file(ssc.str().c_str());
-
- char word[80];
- int mnhc;
- while(file.good())//search and read MNHC
- {
- file>>word;
- if(std::strcmp("MNHC:",word)==0)
- {
- file>>mnhc;
- break;
- }
- else
- {
- file.ignore(150, '\n');
- }
- }
-
- if(mnhc!=MNHC_)
- {
- std::cout<<"please ensure whether 'Restart_md.dat' right!"<>nhcg[i];
- }
- file.get();
- file.ignore(8, '\n');//read eta
- for(int i = 0;i>nhce[i];
- }
- file.get();
- file.ignore(9, '\n');//read peta
- for(int i = 0;i>nhcp[i];
- }
- file.close();
- }
-
- }
-#ifdef __MPI
- MPI_Bcast(&error,1,MPI_INT,0,MPI_COMM_WORLD);
-#endif
- if(error)
- {
- delete[] nhcg;
- delete[] nhcp;
- delete[] nhce;
- exit(0);
- }
-#ifdef __MPI
- MPI_Bcast(nhcg,numIon_*3*MNHC_,MPI_DOUBLE,0,MPI_COMM_WORLD);
- MPI_Bcast(nhcp,numIon_*3*MNHC_,MPI_DOUBLE,0,MPI_COMM_WORLD);
- MPI_Bcast(nhce,numIon_*3*MNHC_,MPI_DOUBLE,0,MPI_COMM_WORLD);
-#endif
-
- for(int i=0;i* vel,
- const double* allmass
- );
- void init_NHC(
- const int &MNHC_in,
- const double &Qmass_in,
- const double &NVT_tau_in,
- const double &dt_in,
- const int &NVT_control,
- std::ofstream &ofs,
- const int &numIon,
- const double &temperature,
- const ModuleBase::Vector3* vel,
- const double* allmass
- );
- double NHChamiltonian(
- const double &KE,
- const double &PE,
- const double &temperature,
- const int &nfrozen
- );
-
- void NHC_info_out(const int& step, const int& recordFreq, const int& numIon);
- void NHC_restart();
-
- private:
- void ADSIntegrator(
- const double &temperature,
- ModuleBase::Vector3* vel,
- const double* allmass
- );
- void LGVIntegrator(
- const double &temperature,
- ModuleBase::Vector3* vel,
- const double* allmass
- );
- void NHCIntegrator(
- const double &temperature,
- ModuleBase::Vector3* vel,
- const double* allmass
- );
-
-
- double gaussrand();
-
- void init_genrand(unsigned long s);
- void init_by_array(unsigned long init_key[], int key_length);
- unsigned long genrand_int32(void);
- long genrand_int31(void);
- double genrand_real1(void);
- double genrand_real2(void);
- double genrand_real3(void);
- double genrand_res53(void);
-
- //NVT thermostat parameters
- const static int N_mt19937=624;
- const static int M_mt19937=397;
- unsigned long const MATRIX_A=0x9908b0dfUL ; /* constant vector a */
- unsigned long const UPPER_MASK=0x80000000UL; /* most significant w-r bits */
- unsigned long const LOWER_MASK=0x7fffffffUL; /* least significant r bits */
- unsigned long mt[N_mt19937]; /* the array for the state vector */
- int mti=625; /* mti==N+1 means mt[N] is not initialized */
- long double gamma; //langevin friction coefficient
- long double nu; //Andersen collision frequency
- long double c1k; //parameter in Langevin
- double c2k; //parameter in Langevin
- const static int nsy=7; //parameter in NHC, constant, no need to modification
- double w[nsy]; //parameters in NHC
- double delta[nsy]; //parameters in NHC
-
- //input parameters
- int MNHC_;
- double Qmass_;
- double NVT_tau_;
- double dt_;
- int numIon_;
-
- //need to be allocated
- ModuleBase::Vector3 *G; //parameter in NHC
- ModuleBase::Vector3 *NHCeta; //NHC position
- ModuleBase::Vector3 *NHCpeta; //NHC momentum
-
-};
-
-#endif
diff --git a/source/module_md/MSST.cpp b/source/module_md/MSST.cpp
new file mode 100644
index 00000000000..d0501608578
--- /dev/null
+++ b/source/module_md/MSST.cpp
@@ -0,0 +1,305 @@
+#include "MSST.h"
+#include "MD_func.h"
+
+MSST::MSST(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para_in, unit_in)
+{
+ std::cout << "MSST" << std::endl;
+
+ mdp.Qmass = mdp.Qmass / pow(ModuleBase::ANGSTROM_AU, 4) / pow(ModuleBase::AU_to_MASS, 2);
+ mdp.velocity = mdp.velocity * ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS;
+ mdp.viscosity = mdp.viscosity / ModuleBase::AU_to_MASS / ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS;
+
+ old_v = new ModuleBase::Vector3 [ucell.nat];
+ dilation.set(1,1,1);
+ omega.set(0,0,0);
+ p0 = 0;
+ e0 = 0;
+ v0 = 1;
+ totmass = 0;
+
+ for(int i=0; i 0 && mdp.tscale > 0)
+ {
+ double fac1 = mdp.tscale * totmass * 2.0 * kinetic / mdp.Qmass;
+ omega[sd] = -1.0 * sqrt(fac1);
+ double fac2 = omega[sd] / v0;
+
+ std::cout << "initial strain rate = " << fac2 << " tscale = " << mdp.tscale << std::endl;
+
+ for(int i=0; i> step_rst_ >> omega[mdp.direction] >> e0 >> v0 >> p0 >> lag_pos;
+
+ file.close();
+ }
+
+#ifdef __MPI
+ MPI_Bcast(&step_rst_, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&omega[mdp.direction], 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&e0, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&v0, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&p0, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&lag_pos, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+#endif
+}
+
+double MSST::extra_term()
+{
+ return 0;
+}
+
+double MSST::vel_sum()
+{
+ double vsum = 0;
+
+ for(int i=0; i const_C = force[i] / allmass[i];
+ ModuleBase::Vector3 const_D;
+ const_D.set(fac/allmass[i], fac/allmass[i], fac/allmass[i]);
+ const_D[sd] -= 2 * omega[sd] / ucell.omega;
+
+ for(int k=0; k<3; ++k)
+ {
+ if( fabs(dthalf*const_D[k]) > 1e-6 )
+ {
+ double expd = exp(dthalf*const_D[k]);
+ vel[i][k] = expd * ( const_C[k] + const_D[k] * vel[i][k] - const_C[k] / expd ) / const_D[k];
+ }
+ else
+ {
+ vel[i][k] += ( const_C[k] + const_D[k] * vel[i][k] ) * dthalf +
+ 0.5 * (const_D[k] * const_D[k] * vel[i][k] + const_C[k] * const_D[k] ) * dthalf * dthalf;
+ }
+ }
+ }
+}
+
+void MSST::propagate_voldot()
+{
+ const int sd = mdp.direction;
+ const double dthalf = 0.5 * mdp.dt;
+ double p_current = stress(sd, sd);
+ double p_msst = mdp.velocity * mdp.velocity * totmass * (v0 - ucell.omega) / (v0 * v0);
+ double const_A = totmass * (p_current - p0 - p_msst) / mdp.Qmass;
+ double const_B = totmass * mdp.viscosity / (mdp.Qmass * ucell.omega);
+
+ // prevent the increase of volume
+ if(ucell.omega > v0 && const_A > 0)
+ {
+ const_A = -const_A;
+ }
+
+ // avoid singularity at B = 0 with Taylor expansion
+ double fac = const_B * dthalf;
+ if(fac > 1e-6)
+ {
+ omega[sd] = (omega[sd] + const_A * (exp(fac) - 1) / const_B) * exp(-fac);
+ }
+ else
+ {
+ omega[sd] += (const_A - const_B * omega[sd]) * dthalf +
+ 0.5 * (const_B * const_B * omega[sd] - const_A * const_B) * dthalf * dthalf;
+ }
+}
\ No newline at end of file
diff --git a/source/module_md/MSST.h b/source/module_md/MSST.h
new file mode 100644
index 00000000000..dfee7a58205
--- /dev/null
+++ b/source/module_md/MSST.h
@@ -0,0 +1,36 @@
+#ifndef MSST_H
+#define MSST_H
+
+#include "verlet.h"
+
+class MSST : public Verlet
+{
+public:
+ MSST(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in);
+ ~MSST();
+
+ void setup();
+ void first_half();
+ void second_half();
+ void outputMD();
+ void write_restart();
+ void restart();
+ double extra_term();
+ double vel_sum();
+ void rescale(double volume);
+ void propagate_vel();
+ void propagate_voldot();
+
+ ModuleBase::Vector3 *old_v;
+ ModuleBase::Vector3 dilation; // dilation scale
+ ModuleBase::Vector3 omega; // time derivative of volume
+ double p0; // initial pressure
+ double v0; // initial volume
+ double e0; // initial energy
+ double totmass; // total mass of the cell
+ double lag_pos; // Lagrangian location of cell
+ double vsum; // sum over v^2
+
+};
+
+#endif
\ No newline at end of file
diff --git a/source/module_md/Makefile b/source/module_md/Makefile
index 698fdb8723e..c9ed7e8fdc8 100644
--- a/source/module_md/Makefile
+++ b/source/module_md/Makefile
@@ -54,7 +54,7 @@ INCLUDES = -I. -Icommands -I${FFTW_INCLUDE_DIR}
#==========================
# -pedantic turns off more extensions and generates more warnings
# -xHost generates instructions for the highest instruction set available on the compilation host processor
-OPTS = ${INCLUDES} -Ofast -std=c++11 -simd -march=native -xHost -m64 -qopenmp -Werror -Wall -pedantic -g
+OPTS = ${INCLUDES} -Ofast -std=c++11 -simd -xHost -m64 -qopenmp -Werror -Wall -pedantic -g
#OPTS = ${INCLUDES} -Ofast -std=c++11 -march=native -Wall -fpermissive -fopenmp
include Makefile.Objects
diff --git a/source/module_md/Makefile.Objects b/source/module_md/Makefile.Objects
index 98c4737105f..36b6e121a38 100644
--- a/source/module_md/Makefile.Objects
+++ b/source/module_md/Makefile.Objects
@@ -7,14 +7,18 @@ HEADERS= *.h
OBJS_MAIN=input.o\
OBJS_MD=driver_classic.o\
-MD_basic.o\
-MD_fire.o\
MD_func.o\
-MD_thermo.o\
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\
OBJS_BASE=matrix3.o\
matrix.o\
diff --git a/source/module_md/NVE.cpp b/source/module_md/NVE.cpp
new file mode 100644
index 00000000000..41e008b6e20
--- /dev/null
+++ b/source/module_md/NVE.cpp
@@ -0,0 +1,51 @@
+#include "NVE.h"
+#include "MD_func.h"
+
+NVE::NVE(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para_in, unit_in){}
+
+NVE::~NVE(){}
+
+void NVE::setup()
+{
+ ModuleBase::TITLE("NVE", "setup");
+ ModuleBase::timer::tick("NVE", "setup");
+
+ Verlet::setup();
+
+ ModuleBase::timer::tick("NVE", "setup");
+}
+
+void NVE::first_half()
+{
+ ModuleBase::TITLE("NVE", "first_half");
+ ModuleBase::timer::tick("NVE", "first_half");
+
+ Verlet::first_half();
+
+ ModuleBase::timer::tick("NVE", "first_half");
+}
+
+void NVE::second_half()
+{
+ ModuleBase::TITLE("NVE", "second_half");
+ ModuleBase::timer::tick("NVE", "second_half");
+
+ Verlet::second_half();
+
+ ModuleBase::timer::tick("NVE", "second_half");
+}
+
+void NVE::outputMD()
+{
+ Verlet::outputMD();
+}
+
+void NVE::write_restart()
+{
+ Verlet::write_restart();
+}
+
+void NVE::restart()
+{
+ Verlet::restart();
+}
\ No newline at end of file
diff --git a/source/module_md/NVE.h b/source/module_md/NVE.h
new file mode 100644
index 00000000000..a71339ffebd
--- /dev/null
+++ b/source/module_md/NVE.h
@@ -0,0 +1,21 @@
+#ifndef NVE_H
+#define NVE_H
+
+#include "verlet.h"
+
+class NVE : public Verlet
+{
+public:
+ NVE(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in);
+ ~NVE();
+
+ void setup();
+ void first_half();
+ void second_half();
+ void outputMD();
+ void write_restart();
+ void restart();
+
+};
+
+#endif
\ No newline at end of file
diff --git a/source/module_md/NVT_ADS.cpp b/source/module_md/NVT_ADS.cpp
new file mode 100644
index 00000000000..54b153f8072
--- /dev/null
+++ b/source/module_md/NVT_ADS.cpp
@@ -0,0 +1,73 @@
+#include "NVT_ADS.h"
+#include "MD_func.h"
+
+NVT_ADS::NVT_ADS(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para_in, unit_in)
+{
+ // convert to a.u. unit
+ mdp.tfreq *= ModuleBase::AU_to_FS;
+
+ nraise = mdp.tfreq * mdp.dt;
+}
+
+NVT_ADS::~NVT_ADS(){}
+
+void NVT_ADS::setup()
+{
+ ModuleBase::TITLE("NVT_ADS", "setup");
+ ModuleBase::timer::tick("NVT_ADS", "setup");
+
+ Verlet::setup();
+
+ ModuleBase::timer::tick("NVT_ADS", "setup");
+}
+
+void NVT_ADS::first_half()
+{
+ ModuleBase::TITLE("NVT_ADS", "first_half");
+ ModuleBase::timer::tick("NVT_ADS", "first_half");
+
+ Verlet::first_half();
+
+ ModuleBase::timer::tick("NVT_ADS", "first_half");
+}
+
+void NVT_ADS::second_half()
+{
+ ModuleBase::TITLE("NVT_ADS", "second_half");
+ ModuleBase::timer::tick("NVT_ADS", "second_half");
+
+ Verlet::second_half();
+
+ double deviation;
+ for(int i=0; i> step_rst_ >> Mnum;
+ if(Mnum != mdp.MNHC)
+ {
+ std::cout<< "Num of NHC is not the same !" << std::endl;
+ ModuleBase::WARNING_QUIT("NVT_NHC", "no Restart_md.dat !");
+ }
+ for(int i=0; i> eta[i];
+ }
+ for(int i=0; i> veta[i];
+ }
+
+ file.close();
+ }
+
+#ifdef __MPI
+ MPI_Bcast(&step_rst_, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(eta, mdp.MNHC, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(veta, mdp.MNHC, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+#endif
+}
+
+void NVT_NHC::integrate()
+{
+ double scale = 1.0;
+ kinetic = MD_func::GetAtomKE(ucell.nat, vel, allmass);
+ double KE = kinetic;
+
+ // update mass
+ update_mass();
+
+ // update force
+ if(Q[0] > 0)
+ {
+ G[0] = (2*KE - (3*ucell.nat - frozen_freedom_)*t_target) / Q[0];
+ }
+ else
+ {
+ G[0] = 0;
+ }
+
+ for(int i=0; i=0; --m)
+ {
+ double aa = exp(-veta[m]*delta/8.0);
+ veta[m] = veta[m] * aa * aa + G[m] * aa * delta /4.0;
+ }
+
+ scale *= exp(-veta[0]*delta/2.0);
+ KE = kinetic * scale * scale;
+
+ // update force
+ if(Q[0] > 0)
+ {
+ G[0] = (2*KE - (3*ucell.nat - frozen_freedom_)*t_target) / Q[0];
+ }
+ else
+ {
+ G[0] = 0;
+ }
+
+ // propogate eta
+ for(int m=0; m CMD_neighbor::cell_periodic(const ModuleBase::Vector
}
//build the neighbor list
-void CMD_neighbor::neighbor(UnitCell_pseudo &ucell_c)
+void CMD_neighbor::neighbor(const UnitCell_pseudo &ucell_c)
{
ModuleBase::TITLE("CMD_neighbor", "Neighbor");
ModuleBase::timer::tick("CMD_neighbor", "Neighbor");
diff --git a/source/module_md/cmd_neighbor.h b/source/module_md/cmd_neighbor.h
index 996892a15dd..af330a3a3de 100644
--- a/source/module_md/cmd_neighbor.h
+++ b/source/module_md/cmd_neighbor.h
@@ -12,7 +12,7 @@ class CMD_neighbor
~CMD_neighbor();
ModuleBase::Vector3 cell_periodic(const ModuleBase::Vector3 a, const ModuleBase::Vector3 b);
- void neighbor(UnitCell_pseudo &ucell_c);
+ void neighbor(const UnitCell_pseudo &ucell_c);
void comm_list(const int num, int *nlist_in, int **list_in, int *nlist_out, int **list_out);
int **list; // record the index of adjent atoms of every atom
diff --git a/source/module_md/driver_classic.cpp b/source/module_md/driver_classic.cpp
index 0b1e2f9b3a9..85812b4f93f 100644
--- a/source/module_md/driver_classic.cpp
+++ b/source/module_md/driver_classic.cpp
@@ -18,10 +18,10 @@ void Driver_classic::init()
ModuleBase::timer::start();
// (1) read the input parameters.
- this->reading();
+ Driver_classic::reading();
// (2) welcome to the classic MD!
- this->classic_world();
+ Driver_classic::classic_world();
// (3) output information
time_t time_finish= std::time(NULL);
@@ -56,6 +56,7 @@ void Driver_classic::convert(UnitCell_pseudo &ucell_c)
if(INPUT.atom_file!="") GlobalV::global_atom_card = INPUT.atom_file;
GlobalV::CALCULATION = INPUT.calculation;
+ GlobalV::BASIS_TYPE = INPUT.basis_type;
GlobalV::OUT_LEVEL = INPUT.out_level;
GlobalV::SEARCH_RADIUS = max(INPUT.search_radius,(INPUT.mdp.rcut_lj+2)*ModuleBase::ANGSTROM_AU);
GlobalV::SEARCH_PBC = INPUT.search_pbc;
@@ -82,7 +83,7 @@ void Driver_classic::classic_world(void)
Run_MD_CLASSIC run_md_classic;
- this->convert(run_md_classic.ucell_c);
+ Driver_classic::convert(run_md_classic.ucell_c);
if(GlobalV::CALCULATION != "md")
{
diff --git a/source/module_md/driver_classic.h b/source/module_md/driver_classic.h
index 4bd9ed21092..b274bae7714 100644
--- a/source/module_md/driver_classic.h
+++ b/source/module_md/driver_classic.h
@@ -10,18 +10,18 @@ class Driver_classic
Driver_classic();
~Driver_classic();
- void init();
+ static void init();
private:
// reading the parameters
- void reading();
+ static void reading();
// convert INPUT parameters for classic MD
- void convert(UnitCell_pseudo &ucell_c);
+ static void convert(UnitCell_pseudo &ucell_c);
// classic MD
- void classic_world();
+ static void classic_world();
};
diff --git a/source/module_md/main.cpp b/source/module_md/main.cpp
index 04195590e7c..ae8708b484c 100644
--- a/source/module_md/main.cpp
+++ b/source/module_md/main.cpp
@@ -9,13 +9,13 @@
int main(int argc, char **argv)
{
+ srand(time(0));
Parallel_Global::read_mpi_parameters(argc,argv);
//----------------------------------------------------------
// main program for doing CMD calculations
//----------------------------------------------------------
- Driver_classic DD;
- DD.init();
+ Driver_classic::init();
#ifdef __MPI
MPI_Finalize();
diff --git a/source/module_md/run_md_classic.cpp b/source/module_md/run_md_classic.cpp
index 4b2927edb40..c13f4831a8a 100644
--- a/source/module_md/run_md_classic.cpp
+++ b/source/module_md/run_md_classic.cpp
@@ -1,31 +1,17 @@
#include "run_md_classic.h"
-#include "MD_basic.h"
-#include "LJ_potential.h"
-#include "DP_potential.h"
-#include "cmd_neighbor.h"
+#include "MD_func.h"
+#include "NVE.h"
+#include "MSST.h"
+#include "FIRE.h"
+#include "NVT_ADS.h"
+#include "NVT_NHC.h"
+#include "Langevin.h"
#include "../input.h"
#include "../src_io/print_info.h"
-#include "../module_neighbor/sltk_atom_arrange.h"
-#include "../module_neighbor/sltk_grid_driver.h"
-Run_MD_CLASSIC::Run_MD_CLASSIC()
-{
- pos_old1 = new double[1];
- pos_old2 = new double[1];
- pos_now = new double[1];
- pos_next = new double[1];
- force = new ModuleBase::Vector3[1];
- stress.create(3,3);
-}
+Run_MD_CLASSIC::Run_MD_CLASSIC(){}
-Run_MD_CLASSIC::~Run_MD_CLASSIC()
-{
- delete[] pos_old1;
- delete[] pos_old2;
- delete[] pos_now;
- delete[] pos_next;
- delete[] force;
-}
+Run_MD_CLASSIC::~Run_MD_CLASSIC(){}
void Run_MD_CLASSIC::classic_md_line(void)
{
@@ -40,162 +26,80 @@ void Run_MD_CLASSIC::classic_md_line(void)
#endif
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
- this->md_allocate_ions();
-
- MD_basic mdb(INPUT.mdp, ucell_c);
- int mdtype = INPUT.mdp.mdtype;
-
- this->istep = 1;
- bool stop = false;
- double potential;
+ // determine the mdtype
+ Verlet *verlet;
+ if(INPUT.mdp.mdtype == -1)
+ {
+ verlet = new FIRE(INPUT.mdp, ucell_c);
+ }
+ else if(INPUT.mdp.mdtype == 0)
+ {
+ verlet = new NVE(INPUT.mdp, ucell_c);
+ }
+ else if(INPUT.mdp.mdtype == 1)
+ {
+ verlet = new NVT_ADS(INPUT.mdp, ucell_c);
+ }
+ else if(INPUT.mdp.mdtype == 2)
+ {
+ verlet = new NVT_NHC(INPUT.mdp, ucell_c);
+ }
+ else if(INPUT.mdp.mdtype == 3)
+ {
+ verlet = new Langevin(INPUT.mdp, ucell_c);
+ }
+ else if(INPUT.mdp.mdtype == 4)
+ {
+ verlet = new MSST(INPUT.mdp, ucell_c);
+ }
- while (istep <= GlobalV::NSTEP && !stop)
+ // md cycle
+ while (verlet->step_ <= GlobalV::NSTEP && !verlet->stop)
{
- time_t fstart = time(NULL);
+ if(verlet->step_ == 0)
+ {
+ verlet->setup();
+ }
+ else
+ {
+ verlet->first_half();
- Print_Info::print_screen(0, 0, istep);
+ // update force and virial due to the update of atom positions
+ MD_func::force_virial(verlet->step_, verlet->mdp, verlet->ucell, verlet->potential, verlet->force, verlet->virial);
- this->md_force_stress(potential);
+ verlet->second_half();
- this->update_pos_classic();
+ MD_func::kinetic_stress(verlet->ucell, verlet->vel, verlet->allmass, verlet->kinetic, verlet->stress);
- if (mdtype == 1 || mdtype == 2)
- {
- mdb.runNVT(istep, potential, force, stress);
- }
- else if (mdtype == 0)
- {
- mdb.runNVE(istep, potential, force, stress);
- }
- else if (mdtype == -1)
- {
- stop = mdb.runFIRE(istep, potential, force, stress);
+ verlet->stress += verlet->virial;
}
- else
+
+ if(verlet->step_ % verlet->mdp.recordFreq == 0)
{
- ModuleBase::WARNING_QUIT("md_cells_classic", "mdtype should be -1~2!");
+ Print_Info::print_screen(0, 0, verlet->step_ + verlet->step_rst_);
+ verlet->outputMD();
+
+ verlet->ucell.update_vel(verlet->vel);
+ std::stringstream file;
+ file << GlobalV::global_out_dir << "STRU_MD_" << verlet->step_ + verlet->step_rst_;
+#ifdef __LCAO
+ verlet->ucell.print_stru_file(GlobalC::ORB, file.str(), 1, 1);
+#else
+ verlet->ucell.print_stru_file(file.str(), 1, 1);
+#endif
+ MD_func::MDdump(verlet->step_ + verlet->step_rst_, verlet->ucell.nat, verlet->virial, verlet->force);
+ verlet->write_restart();
}
- time_t fend = time(NULL);
-
- ucell_c.save_cartesian_position(this->pos_next);
-
- ++istep;
+ verlet->step_++;
}
GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
GlobalV::ofs_running << std::setprecision(16);
- GlobalV::ofs_running << " !FINAL_ETOT_IS " << potential*ModuleBase::Hartree_to_eV << " eV" << std::endl;
+ GlobalV::ofs_running << " !FINAL_ETOT_IS " << verlet->potential*ModuleBase::Hartree_to_eV << " eV" << std::endl;
GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
- ModuleBase::timer::tick("Run_MD_CLASSIC", "md_cells_classic");
+ delete verlet;
+ ModuleBase::timer::tick("Run_MD_CLASSIC", "classic_md_line");
return;
-}
-
-void Run_MD_CLASSIC::md_force_stress(double &potential)
-{
- ModuleBase::TITLE("Run_MD_CLASSIC", "md_force_stress");
- ModuleBase::timer::tick("Run_MD_CLASSIC", "md_force_stress");
-
- if(INPUT.mdp.md_potential == "LJ")
- {
- bool which_method = this->ucell_c.judge_big_cell();
- if(which_method)
- {
- CMD_neighbor cmd_neigh;
- cmd_neigh.neighbor(this->ucell_c);
-
- potential = LJ_potential::Lennard_Jones(
- this->ucell_c,
- cmd_neigh,
- this->force,
- this->stress);
- }
- else
- {
- Grid_Driver grid_neigh(GlobalV::test_deconstructor, GlobalV::test_grid_driver, GlobalV::test_grid);
- atom_arrange::search(
- GlobalV::SEARCH_PBC,
- GlobalV::ofs_running,
- grid_neigh,
- this->ucell_c,
- GlobalV::SEARCH_RADIUS,
- GlobalV::test_atom_input,
- INPUT.test_just_neighbor);
-
- potential = LJ_potential::Lennard_Jones(
- this->ucell_c,
- grid_neigh,
- this->force,
- this->stress);
- }
- }
- else if(INPUT.mdp.md_potential == "DP")
- {
- DP_potential::DP_pot(ucell_c, potential, force, stress);
- }
- else if(INPUT.mdp.md_potential == "FP")
- {
- ModuleBase::WARNING_QUIT("md_force_stress", "FPMD is only available in integrated program or PW module !");
- }
- else
- {
- ModuleBase::WARNING_QUIT("md_force_stress", "Unsupported MD potential !");
- }
-
- ModuleBase::GlobalFunc::NEW_PART(" TOTAL-FORCE (eV/Angstrom)");
- GlobalV::ofs_running << std::endl;
- GlobalV::ofs_running << " atom x y z" << std::endl;
- const double fac = ModuleBase::Hartree_to_eV*ModuleBase::ANGSTROM_AU;
- int iat = 0;
- for(int it=0; itpos_old1;
- delete[] this->pos_old2;
- delete[] this->pos_now;
- delete[] this->pos_next;
- delete[] this->force;
-
- 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];
- this->force = new ModuleBase::Vector3[ucell_c.nat];
-
- ModuleBase::GlobalFunc::ZEROS(pos_old1, pos_dim);
- ModuleBase::GlobalFunc::ZEROS(pos_old2, pos_dim);
- ModuleBase::GlobalFunc::ZEROS(pos_now, pos_dim);
- ModuleBase::GlobalFunc::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];
- }
- this->ucell_c.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 0f79ddfabb4..cd675c112cd 100644
--- a/source/module_md/run_md_classic.h
+++ b/source/module_md/run_md_classic.h
@@ -14,20 +14,7 @@ class Run_MD_CLASSIC
void classic_md_line(void);
- void md_force_stress(double &potential);
- void md_allocate_ions(void);
- void update_pos_classic(void);
-
- ModuleBase::Vector3 *force; //force of each atom
- ModuleBase::matrix stress; //stress for this lattice
-
-private:
- int istep;
- double* pos_old1;
- double* pos_old2;
- double* pos_now;
- double* pos_next;
- int pos_dim;
+
};
#endif
\ No newline at end of file
diff --git a/source/module_md/verlet.cpp b/source/module_md/verlet.cpp
new file mode 100644
index 00000000000..594418f754a
--- /dev/null
+++ b/source/module_md/verlet.cpp
@@ -0,0 +1,185 @@
+#include "verlet.h"
+#include "MD_func.h"
+
+Verlet::Verlet(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in):
+ mdp(MD_para_in),
+ ucell(unit_in)
+{
+ stop = false;
+
+ allmass = new double [ucell.nat];
+ pos = new ModuleBase::Vector3 [ucell.nat];
+ vel = new ModuleBase::Vector3 [ucell.nat];
+ ionmbl = new ModuleBase::Vector3 [ucell.nat];
+ force = new ModuleBase::Vector3 [ucell.nat];
+ virial.create(3,3);
+ stress.create(3,3);
+
+ // convert to a.u. unit
+ mdp.dt /= ModuleBase::AU_to_FS;
+ mdp.tfirst /= ModuleBase::Hartree_to_K;
+ mdp.tlast /= ModuleBase::Hartree_to_K;
+
+ // LJ parameters
+ mdp.rcut_lj *= ModuleBase::ANGSTROM_AU;
+ mdp.epsilon_lj /= ModuleBase::Hartree_to_eV;
+ mdp.sigma_lj *= ModuleBase::ANGSTROM_AU;
+
+ step_ = 0;
+ step_rst_ = 0;
+ if(mdp.rstMD) unit_in.set_vel = 1;
+
+ MD_func::InitPos(ucell, pos);
+ MD_func::InitVel(ucell, mdp.tfirst, allmass, frozen_freedom_, ionmbl, vel);
+}
+
+Verlet::~Verlet()
+{
+ delete []allmass;
+ delete []pos;
+ delete []vel;
+ delete []ionmbl;
+ delete []force;
+}
+
+void Verlet::setup()
+{
+ if(mdp.rstMD)
+ {
+ restart();
+ }
+
+ MD_func::force_virial(step_, mdp, ucell, potential, force, virial);
+ MD_func::kinetic_stress(ucell, vel, allmass, kinetic, stress);
+ stress += virial;
+
+ temperature_ = 2*kinetic/(double(3*ucell.nat-frozen_freedom_))*ModuleBase::Hartree_to_K;
+}
+
+void Verlet::first_half()
+{
+ for(int i=0; i> step_rst_;
+
+ file.close();
+ }
+
+#ifdef __MPI
+ MPI_Bcast(&step_rst_, 1, MPI_INT, 0, MPI_COMM_WORLD);
+#endif
+}
\ No newline at end of file
diff --git a/source/module_md/verlet.h b/source/module_md/verlet.h
new file mode 100644
index 00000000000..d07910e5d87
--- /dev/null
+++ b/source/module_md/verlet.h
@@ -0,0 +1,43 @@
+#ifndef VERLET_H
+#define VERLET_H
+
+#include "MD_parameters.h"
+#include "../module_cell/unitcell_pseudo.h"
+
+class Verlet
+{
+public:
+ Verlet(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in);
+ virtual ~Verlet();
+
+ virtual void setup();
+ virtual void first_half();
+ virtual void second_half();
+ virtual void outputMD();
+ virtual void write_restart();
+ virtual void restart();
+
+ MD_parameters &mdp;
+ UnitCell_pseudo &ucell;
+ bool stop; // MD stop or not
+
+ // All parameters are in a.u. unit.
+ double temperature_;
+ int step_;
+ int step_rst_;
+ double energy_;
+ int frozen_freedom_;
+
+ double *allmass; // atom mass
+ ModuleBase::Vector3 *pos; // atom position
+ ModuleBase::Vector3 *vel; // atom velocity
+ ModuleBase::Vector3 *ionmbl; // atom is frozen or not
+ ModuleBase::Vector3 *force; // force of each atom
+ ModuleBase::matrix virial; // virial for this lattice
+ ModuleBase::matrix stress; // stress for this lattice
+ double potential; // potential energy
+ double kinetic; // kinetic energy
+
+};
+
+#endif
\ No newline at end of file
diff --git a/source/module_orbital/ORB_read.cpp b/source/module_orbital/ORB_read.cpp
index 028fb93f999..5f5a8ee4128 100644
--- a/source/module_orbital/ORB_read.cpp
+++ b/source/module_orbital/ORB_read.cpp
@@ -556,7 +556,7 @@ void LCAO_Orbitals::read_orb_file(
this->kmesh,
this->dk,
this->dr_uniform,
- true,
+ GlobalV::out_element_info,
true,
force_flag); // delta k mesh in reciprocal space
diff --git a/source/module_symmetry/symmetry.cpp b/source/module_symmetry/symmetry.cpp
index 8df9de04e5c..24524428541 100644
--- a/source/module_symmetry/symmetry.cpp
+++ b/source/module_symmetry/symmetry.cpp
@@ -20,7 +20,7 @@ Symmetry::~Symmetry()
bool Symmetry::symm_flag=false;
-void Symmetry::analy_sys(const UnitCell_pseudo &ucell, const output &out, std::ofstream &ofs_running)
+void Symmetry::analy_sys(const UnitCell_pseudo &ucell, std::ofstream &ofs_running)
{
if (available == false) return;
ModuleBase::TITLE("Symmetry","init");
@@ -75,7 +75,7 @@ void Symmetry::analy_sys(const UnitCell_pseudo &ucell, const output &out, std::o
// std::cout << "a1 = " << a2.x << " " << a2.y << " " << a2.z < s1, s2, s3;
diff --git a/source/run_lcao.cpp b/source/run_lcao.cpp
index 1c118270781..965898bf7a3 100644
--- a/source/run_lcao.cpp
+++ b/source/run_lcao.cpp
@@ -23,9 +23,9 @@ void Run_lcao::lcao_line(void)
// improvement: a) separating the first reading of the atom_card and subsequent
// cell relaxation. b) put GlobalV::NLOCAL and GlobalV::NBANDS as input parameters
#ifdef __LCAO
- GlobalC::ucell.setup_cell( GlobalC::ORB, GlobalV::global_pseudo_dir, GlobalC::out, GlobalV::global_atom_card, GlobalV::ofs_running);
+ GlobalC::ucell.setup_cell( GlobalC::ORB, GlobalV::global_pseudo_dir, GlobalV::global_atom_card, GlobalV::ofs_running);
#else
- GlobalC::ucell.setup_cell( GlobalV::global_pseudo_dir, GlobalC::out, GlobalV::global_atom_card, GlobalV::ofs_running);
+ GlobalC::ucell.setup_cell( GlobalV::global_pseudo_dir, GlobalV::global_atom_card, GlobalV::ofs_running);
#endif
if(INPUT.test_just_neighbor)
{
@@ -66,7 +66,7 @@ void Run_lcao::lcao_line(void)
// symmetry analysis should be performed every time the cell is changed
if (ModuleSymmetry::Symmetry::symm_flag)
{
- GlobalC::symm.analy_sys(GlobalC::ucell, GlobalC::out, GlobalV::ofs_running);
+ GlobalC::symm.analy_sys(GlobalC::ucell, GlobalV::ofs_running);
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
}
diff --git a/source/run_pw.cpp b/source/run_pw.cpp
index c94b4a7f105..4efc6ed33fe 100644
--- a/source/run_pw.cpp
+++ b/source/run_pw.cpp
@@ -24,9 +24,9 @@ void Run_pw::plane_wave_line(void)
// improvement: a) separating the first reading of the atom_card and subsequent
// cell relaxation. b) put GlobalV::NLOCAL and GlobalV::NBANDS as input parameters
#ifdef __LCAO
- GlobalC::ucell.setup_cell( GlobalC::ORB, GlobalV::global_pseudo_dir, GlobalC::out, GlobalV::global_atom_card, GlobalV::ofs_running);
+ GlobalC::ucell.setup_cell( GlobalC::ORB, GlobalV::global_pseudo_dir, GlobalV::global_atom_card, GlobalV::ofs_running);
#else
- GlobalC::ucell.setup_cell( GlobalV::global_pseudo_dir, GlobalC::out, GlobalV::global_atom_card, GlobalV::ofs_running);
+ GlobalC::ucell.setup_cell( GlobalV::global_pseudo_dir, GlobalV::global_atom_card, GlobalV::ofs_running);
#endif
//GlobalC::ucell.setup_cell( GlobalV::global_pseudo_dir , GlobalV::global_atom_card , GlobalV::ofs_running, GlobalV::NLOCAL, GlobalV::NBANDS);
@@ -49,7 +49,7 @@ void Run_pw::plane_wave_line(void)
// symmetry analysis should be performed every time the cell is changed
if (ModuleSymmetry::Symmetry::symm_flag)
{
- GlobalC::symm.analy_sys(GlobalC::ucell, GlobalC::out, GlobalV::ofs_running);
+ GlobalC::symm.analy_sys(GlobalC::ucell, GlobalV::ofs_running);
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
}
diff --git a/source/src_io/CMakeLists.txt b/source/src_io/CMakeLists.txt
index aeef1109c7f..f7dc7e1a70b 100644
--- a/source/src_io/CMakeLists.txt
+++ b/source/src_io/CMakeLists.txt
@@ -39,4 +39,5 @@ add_library(
write_rho.cpp
write_rho_cube.cpp
write_rho_dipole.cpp
+ write_wfc_realspace.cpp
)
diff --git a/source/src_io/dos.cpp b/source/src_io/dos.cpp
index 773f3a6bf60..2169dba0a38 100644
--- a/source/src_io/dos.cpp
+++ b/source/src_io/dos.cpp
@@ -141,14 +141,14 @@ void Dos::calculate_Mulliken(const std::string &fa)
bool Dos::calculate_dos
(
const int &is,
- const int *isk,
+ const std::vector &isk,
const std::string &fa, //file address
const double &de_ev, // delta energy in ev
const double &emax_ev,
const double &emin_ev,// minimal energy in ev.
const int &nks,//number of k points
const int &nkstot,
- const double *wk,//weight of k points
+ const std::vector &wk,//weight of k points
const ModuleBase::matrix &wg,//weight of (kpoint,bands)
const int &nbands,// number of bands
double** ekb//store energy for each k point and each band
diff --git a/source/src_io/dos.h b/source/src_io/dos.h
index 9662b737404..1194e12b1c2 100644
--- a/source/src_io/dos.h
+++ b/source/src_io/dos.h
@@ -6,14 +6,14 @@ namespace Dos
{
bool calculate_dos(
const int &is,
- const int *isk,
+ const std::vector &isk,
const std::string &fn,// file address.
const double &de_ev, // delta energy in ev.
const double &emax_ev,// maximal energy in ev.
const double &emin_ev,// minimal energy in ev.
const int &nks,//number of k points
const int &nkstot,
- const double *wk,//weight of k points
+ const std::vector &wk,//weight of k points
const ModuleBase::matrix &wg,//weight of (kpoint,bands)
const int &nbands,// number of bands
double **ekb);//store energy for each k point and each band
diff --git a/source/src_io/numerical_basis.cpp b/source/src_io/numerical_basis.cpp
index 5525e14b18d..5013ec693d2 100644
--- a/source/src_io/numerical_basis.cpp
+++ b/source/src_io/numerical_basis.cpp
@@ -67,7 +67,7 @@ void Numerical_Basis::output_overlap( const ModuleBase::ComplexMatrix *psi)
std::ofstream ofs;
std::stringstream ss;
// the parameter 'winput::spillage_outdir' is read from INPUTw.
- ss << winput::spillage_outdir << "/" << GlobalC::ucell.latName << "." << derivative_order << ".dat";
+ ss << winput::spillage_outdir << "/" << "orb_matrix." << derivative_order << ".dat";
if (GlobalV::MY_RANK==0)
{
ofs.open(ss.str().c_str());
diff --git a/source/src_io/output.cpp b/source/src_io/output.cpp
index 2df23801f9c..defe5b1cebe 100644
--- a/source/src_io/output.cpp
+++ b/source/src_io/output.cpp
@@ -2,7 +2,7 @@
#include "output.h"
-void output::printrm(std::ofstream &ofs,const std::string &s, const ModuleBase::matrix &m, const double &limit) const
+void output::printrm(std::ofstream &ofs,const std::string &s, const ModuleBase::matrix &m, const double &limit)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -29,7 +29,7 @@ void output::printrm(std::ofstream &ofs,const std::string &s, const ModuleBase::
return;
}
-void output::printrm(const std::string &s, const ModuleBase::matrix &m, const double &limit) const
+void output::printrm(const std::string &s, const ModuleBase::matrix &m, const double &limit)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -55,7 +55,7 @@ void output::printrm(const std::string &s, const ModuleBase::matrix &m, const do
return;
}
-void output::printcm_norm(std::ofstream &ofs, const std::string &s, const ModuleBase::ComplexMatrix &m, const double &limit)const
+void output::printcm_norm(std::ofstream &ofs, const std::string &s, const ModuleBase::ComplexMatrix &m, const double &limit)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -75,7 +75,7 @@ void output::printcm_norm(std::ofstream &ofs, const std::string &s, const Module
return;
}
-void output::printcm_norm(const std::string &s, const ModuleBase::ComplexMatrix &m, const double &limit)const
+void output::printcm_norm(const std::string &s, const ModuleBase::ComplexMatrix &m, const double &limit)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -96,7 +96,7 @@ void output::printcm_norm(const std::string &s, const ModuleBase::ComplexMatrix
}
-void output::printcm(std::ofstream &ofs, const std::string &s, const ModuleBase::ComplexMatrix &m) const
+void output::printcm(std::ofstream &ofs, const std::string &s, const ModuleBase::ComplexMatrix &m)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -117,7 +117,7 @@ void output::printcm(std::ofstream &ofs, const std::string &s, const ModuleBase:
return;
}//end print cm
-void output::printcm(const std::string &s, const ModuleBase::ComplexMatrix &m) const
+void output::printcm(const std::string &s, const ModuleBase::ComplexMatrix &m)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -136,7 +136,7 @@ void output::printcm(const std::string &s, const ModuleBase::ComplexMatrix &m) c
return;
}
-void output::printcm_real_limit_hermit(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit) const
+void output::printcm_real_limit_hermit(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -157,7 +157,7 @@ void output::printcm_real_limit_hermit(const std::string &s, const ModuleBase::C
return;
}
-void output::printcm_real(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit) const
+void output::printcm_real(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -180,7 +180,7 @@ void output::printcm_real(const std::string &s, const ModuleBase::ComplexMatrix
-void output::printcm_imag(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit) const
+void output::printcm_imag(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit)
{
const int b1 = m.nr;
const int b2 = m.nc;
@@ -205,7 +205,7 @@ void output::printcm_imag(const std::string &s, const ModuleBase::ComplexMatrix
-void output::printr3_d(std::ofstream &ofs, const std::string &s,const ModuleBase::realArray &u) const
+void output::printr3_d(std::ofstream &ofs, const std::string &s,const ModuleBase::realArray &u)
{
const int b1 = u.getBound1();
const int b2 = u.getBound2();
@@ -228,7 +228,7 @@ void output::printr3_d(std::ofstream &ofs, const std::string &s,const ModuleBase
return;
}//end printr3_d
-void output::printr4_d(std::ofstream &ofs, const std::string &s,const ModuleBase::realArray &u) const
+void output::printr4_d(std::ofstream &ofs, const std::string &s,const ModuleBase::realArray &u)
{
const int b1 = u.getBound1();
const int b2 = u.getBound2();
@@ -252,7 +252,7 @@ void output::printr4_d(std::ofstream &ofs, const std::string &s,const ModuleBase
}
}//end print4_d
-void output::printM3(std::ofstream &ofs,const std::string &description, const ModuleBase::Matrix3 &m)const
+void output::printM3(std::ofstream &ofs,const std::string &description, const ModuleBase::Matrix3 &m)
{
ofs << " " << description << std::endl;
ofs << std::setiosflags(std::ios::showpos);
@@ -263,7 +263,7 @@ void output::printM3(std::ofstream &ofs,const std::string &description, const Mo
return;
}
-void output::printM3(const std::string &description, const ModuleBase::Matrix3 &m)const
+void output::printM3(const std::string &description, const ModuleBase::Matrix3 &m)
{
std::cout << "\n " << description << std::endl;
std::cout << std::setw(20) << m.e11 << std::setw(20) << m.e12 << std::setw(20) << m.e13 << "\n"
diff --git a/source/src_io/output.h b/source/src_io/output.h
index 752931c9052..3a4663c83d6 100644
--- a/source/src_io/output.h
+++ b/source/src_io/output.h
@@ -15,36 +15,36 @@ class output
//============================
// Print realArray (3D or 4D)
//============================
- void printr3_d(std::ofstream &ofs,const std::string &s,const ModuleBase::realArray &u) const;
- void printr4_d(std::ofstream &ofs,const std::string &s,const ModuleBase::realArray &u) const;
+ static void printr3_d(std::ofstream &ofs,const std::string &s,const ModuleBase::realArray &u);
+ static void printr4_d(std::ofstream &ofs,const std::string &s,const ModuleBase::realArray &u);
//===========================
// print matrix3
//===========================
- void printM3(std::ofstream &ofs,const std::string& description, const ModuleBase::Matrix3 &m)const;
- void printM3(const std::string &description, const ModuleBase::Matrix3 &m)const;
+ static void printM3(std::ofstream &ofs,const std::string& description, const ModuleBase::Matrix3 &m);
+ static void printM3(const std::string &description, const ModuleBase::Matrix3 &m);
//===============================
// print matrix
//===============================
- void printrm(std::ofstream &ofs,const std::string &s, const ModuleBase::matrix &m, const double &limit = 1.0e-15) const;
- void printrm(const std::string &s, const ModuleBase::matrix &m, const double &limit = 1.0e-15) const;
+ static void printrm(std::ofstream &ofs,const std::string &s, const ModuleBase::matrix &m, const double &limit = 1.0e-15);
+ static void printrm(const std::string &s, const ModuleBase::matrix &m, const double &limit = 1.0e-15);
//===============================
// print ModuleBase::ComplexMatrix
//===============================
- void printcm(std::ofstream &ofs,const std::string &s, const ModuleBase::ComplexMatrix &m) const;
+ static void printcm(std::ofstream &ofs,const std::string &s, const ModuleBase::ComplexMatrix &m);
- void printcm(const std::string &s, const ModuleBase::ComplexMatrix &m) const;
+ static void printcm(const std::string &s, const ModuleBase::ComplexMatrix &m);
- void printcm_real(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit = 1.0e-15) const;
+ static void printcm_real(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit = 1.0e-15);
- void printcm_real_limit_hermit(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit) const;
+ static void printcm_real_limit_hermit(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit);
- void printcm_imag(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit = 1.0e-15) const;
- void printcm_norm(const std::string &s, const ModuleBase::ComplexMatrix &m, const double &limit)const;
- void printcm_norm(std::ofstream &ofs, const std::string &s, const ModuleBase::ComplexMatrix &m, const double &limit)const;
+ static void printcm_imag(const std::string &s, const ModuleBase::ComplexMatrix &m,const double &limit = 1.0e-15);
+ static void printcm_norm(const std::string &s, const ModuleBase::ComplexMatrix &m, const double &limit);
+ static void printcm_norm(std::ofstream &ofs, const std::string &s, const ModuleBase::ComplexMatrix &m, const double &limit);
//***************
@@ -53,7 +53,7 @@ class output
public:
template
- void printr1_d(std::ofstream &ofs, const std::string &s,T *u, int n1) const
+ static void printr1_d(std::ofstream &ofs, const std::string &s,T *u, int n1)
{
ofs<<"\n\n "<
- void printV3(const ModuleBase::Vector3 v)const
+ static void printV3(const ModuleBase::Vector3 v)
{
std::cout << " ";
std::cout << std::setw(18) << v.x << std::setw(18) << v.y << std::setw(18) << v.z << std::endl;
}
template
- void printv31_d(std::ofstream &ofs, const std::string &s, ModuleBase::Vector3 *u, int n1) const
+ static void printv31_d(std::ofstream &ofs, const std::string &s, ModuleBase::Vector3 *u, int n1)
{
ofs << " " << s << " Dimension = " << n1 << std::endl;
if (n1 <= 0)return;
@@ -117,7 +117,7 @@ class output
}
template
- void printv31_d(const std::string &s, ModuleBase::Vector3 *u, int n1) const
+ static void printv31_d(const std::string &s, ModuleBase::Vector3 *u, int n1)
{
std::cout << "\n " << s << " dimension = " << n1;
if (n1 <= 0)return;
diff --git a/source/src_io/print_info.cpp b/source/src_io/print_info.cpp
index 453ede338fa..da722bf9a0e 100644
--- a/source/src_io/print_info.cpp
+++ b/source/src_io/print_info.cpp
@@ -40,7 +40,6 @@ void Print_Info::setup_parameters(UnitCell_pseudo &ucell, K_Vectors &kv, xcfunc
if(INPUT.mdp.mdtype ==1 || INPUT.mdp.mdtype==2)
{
std::cout << " ENSEMBLE : " << "NVT" << std::endl;
- std::cout << " Qmass for NVT(a.u.) : " << INPUT.mdp.Qmass/ModuleBase::AU_to_MASS << std::endl;
}
else if(INPUT.mdp.mdtype==0)
{
diff --git a/source/src_io/write_HS.cpp b/source/src_io/write_HS.cpp
index f97846d3466..73950357153 100644
--- a/source/src_io/write_HS.cpp
+++ b/source/src_io/write_HS.cpp
@@ -46,7 +46,7 @@ void HS_Matrix::save_HS_ccf(const int &iter, const int &Hnnz, const int *colptr_
}
else
{
- // mohan update 2021-02-10
+ // mohan update 2021-02-10
ssh << GlobalV::global_out_dir << "H" << ELEC_scf::iter << "_" << iter+1 << ".ccf";
sss << GlobalV::global_out_dir << "S" << ELEC_scf::iter << "_" << iter+1 << ".ccf";
}
@@ -861,94 +861,138 @@ void HS_Matrix::save_HSR_tr(const int current_spin)
return;
}
-void HS_Matrix::save_HSR_sparse(const int ¤t_spin, const double &sparse_threshold, const bool &binary)
+void HS_Matrix::save_HSR_sparse(
+ const double &sparse_threshold,
+ const bool &binary,
+ const std::string &SR_filename,
+ const std::string &HR_filename_up,
+ const std::string &HR_filename_down = ""
+)
{
ModuleBase::TITLE("HS_Matrix","save_HSR_sparse");
ModuleBase::timer::tick("HS_Matrix","save_HSR_sparse");
+ auto &all_R_coor_ptr = GlobalC::LM.all_R_coor;
auto &HR_sparse_ptr = GlobalC::LM.HR_sparse;
- auto &HR_soc_sparse_ptr = GlobalC::LM.HR_soc_sparse;
auto &SR_sparse_ptr = GlobalC::LM.SR_sparse;
+ auto &HR_soc_sparse_ptr = GlobalC::LM.HR_soc_sparse;
auto &SR_soc_sparse_ptr = GlobalC::LM.SR_soc_sparse;
- int R_x = GlobalC::GridD.getCellX();
- int R_y = GlobalC::GridD.getCellY();
- int R_z = GlobalC::GridD.getCellZ();
+ int total_R_num = all_R_coor_ptr.size();
+ int output_R_number = 0;
+ int *H_nonzero_num[2] = {nullptr, nullptr};
+ int *S_nonzero_num = nullptr;
- double R_minX = GlobalC::GridD.getD_minX();
- double R_minY = GlobalC::GridD.getD_minY();
- double R_minZ = GlobalC::GridD.getD_minZ();
+ S_nonzero_num = new int[total_R_num];
+ ModuleBase::GlobalFunc::ZEROS(S_nonzero_num, total_R_num);
- int total_R_number = R_x * R_y * R_z;
- int output_R_number = 0;
- int *H_nonzero_number = new int[total_R_number];
- int *S_nonzero_number = new int[total_R_number];
- int count_n = 0;
- for (int ix = 0; ix < R_x; ++ix)
+ int spin_loop = 1;
+ if (GlobalV::NSPIN == 2)
+ {
+ spin_loop = 2;
+ }
+
+ for (int ispin = 0; ispin < spin_loop; ++ispin)
{
- for (int iy = 0; iy < R_y; ++iy)
+ H_nonzero_num[ispin] = new int[total_R_num];
+ ModuleBase::GlobalFunc::ZEROS(H_nonzero_num[ispin], total_R_num);
+ }
+
+ int count = 0;
+ for (auto &R_coor : all_R_coor_ptr)
+ {
+ if (GlobalV::NSPIN != 4)
{
- for (int iz = 0; iz < R_z; ++iz)
+ for (int ispin = 0; ispin < spin_loop; ++ispin)
{
- H_nonzero_number[count_n] = 0;
- S_nonzero_number[count_n] = 0;
- if (GlobalV::NSPIN != 4)
+ auto iter = HR_sparse_ptr[ispin].find(R_coor);
+ if (iter != HR_sparse_ptr[ispin].end())
{
- for (auto &iter : HR_sparse_ptr[ix][iy][iz])
- {
- H_nonzero_number[count_n] += iter.second.size();
- }
- for (auto &iter : SR_sparse_ptr[ix][iy][iz])
+ for (auto &row_loop : iter->second)
{
- S_nonzero_number[count_n] += iter.second.size();
+ H_nonzero_num[ispin][count] += row_loop.second.size();
}
}
- else
+ }
+
+ auto iter = SR_sparse_ptr.find(R_coor);
+ if (iter != SR_sparse_ptr.end())
+ {
+ for (auto &row_loop : iter->second)
{
- for (auto &iter : HR_soc_sparse_ptr[ix][iy][iz])
- {
- H_nonzero_number[count_n] += iter.second.size();
- }
- for (auto &iter : SR_soc_sparse_ptr[ix][iy][iz])
- {
- S_nonzero_number[count_n] += iter.second.size();
- }
+ S_nonzero_num[count] += row_loop.second.size();
}
+ }
+ }
+ else
+ {
+ auto iter = HR_soc_sparse_ptr.find(R_coor);
+ if (iter != HR_soc_sparse_ptr.end())
+ {
+ for (auto &row_loop : iter->second)
+ {
+ H_nonzero_num[0][count] += row_loop.second.size();
+ }
+ }
- count_n++;
+ iter = SR_soc_sparse_ptr.find(R_coor);
+ if (iter != SR_soc_sparse_ptr.end())
+ {
+ for (auto &row_loop : iter->second)
+ {
+ S_nonzero_num[count] += row_loop.second.size();
+ }
}
}
+
+ count++;
}
- Parallel_Reduce::reduce_int_all(H_nonzero_number, total_R_number);
- Parallel_Reduce::reduce_int_all(S_nonzero_number, total_R_number);
+ Parallel_Reduce::reduce_int_all(S_nonzero_num, total_R_num);
+ for (int ispin = 0; ispin < spin_loop; ++ispin)
+ {
+ Parallel_Reduce::reduce_int_all(H_nonzero_num[ispin], total_R_num);
+ }
- for (int index = 0; index < total_R_number; ++index)
+ if (GlobalV::NSPIN == 2)
{
- if (H_nonzero_number[index] == 0 && S_nonzero_number[index] == 0)
+ for (int index = 0; index < total_R_num; ++index)
{
- // do nothing
- }
- else
+ if (H_nonzero_num[0][index] != 0 || H_nonzero_num[1][index] != 0 || S_nonzero_num[index] != 0)
+ {
+ output_R_number++;
+ }
+ }
+ }
+ else
+ {
+ for (int index = 0; index < total_R_num; ++index)
{
- output_R_number++;
+ if (H_nonzero_num[0][index] != 0 || S_nonzero_num[index] != 0)
+ {
+ output_R_number++;
+ }
}
}
- std::stringstream ssh;
+ std::stringstream ssh[2];
std::stringstream sss;
- ssh << GlobalV::global_out_dir << "data-HR-sparse_SPIN" << current_spin << ".csr";
- sss << GlobalV::global_out_dir << "data-SR-sparse_SPIN" << current_spin << ".csr";
- std::ofstream g1;
+ ssh[0] << GlobalV::global_out_dir << HR_filename_up;
+ ssh[1] << GlobalV::global_out_dir << HR_filename_down;
+ sss << GlobalV::global_out_dir << SR_filename;
+ std::ofstream g1[2];
std::ofstream g2;
if(GlobalV::DRANK==0)
{
if (binary)
{
- g1.open(ssh.str().c_str(), ios::binary);
- g1.write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int));
- g1.write(reinterpret_cast(&output_R_number), sizeof(int));
+ for (int ispin = 0; ispin < spin_loop; ++ispin)
+ {
+ g1[ispin].open(ssh[ispin].str().c_str(), ios::binary);
+ g1[ispin].write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int));
+ g1[ispin].write(reinterpret_cast(&output_R_number), sizeof(int));
+ }
g2.open(sss.str().c_str(), ios::binary);
g2.write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int));
@@ -956,9 +1000,12 @@ void HS_Matrix::save_HSR_sparse(const int ¤t_spin, const double &sparse_th
}
else
{
- g1.open(ssh.str().c_str());
- g1 << "Matrix Dimension of H(R): " << GlobalV::NLOCAL <(&dRx), sizeof(int));
+ g1[ispin].write(reinterpret_cast(&dRy), sizeof(int));
+ g1[ispin].write(reinterpret_cast(&dRz), sizeof(int));
+ g1[ispin].write(reinterpret_cast(&H_nonzero_num[ispin][count]), sizeof(int));
}
- if (GlobalV::DRANK == 0)
+ g2.write(reinterpret_cast(&dRx), sizeof(int));
+ g2.write(reinterpret_cast(&dRy), sizeof(int));
+ g2.write(reinterpret_cast(&dRz), sizeof(int));
+ g2.write(reinterpret_cast(&S_nonzero_num[count]), sizeof(int));
+ }
+ else
+ {
+ for (int ispin = 0; ispin < spin_loop; ++ispin)
{
- if (binary)
- {
- g1.write(reinterpret_cast(&dRx), sizeof(int));
- g1.write(reinterpret_cast(&dRy), sizeof(int));
- g1.write(reinterpret_cast(&dRz), sizeof(int));
- g1.write(reinterpret_cast(&H_nonzero_number[count_n]), sizeof(int));
-
- g2.write(reinterpret_cast(&dRx), sizeof(int));
- g2.write(reinterpret_cast(&dRy), sizeof(int));
- g2.write(reinterpret_cast(&dRz), sizeof(int));
- g2.write(reinterpret_cast(&S_nonzero_number[count_n]), sizeof(int));
- }
- else
- {
- g1 << dRx << " " << dRy << " " << dRz << " " << H_nonzero_number[count_n] << std::endl;
- g2 << dRx << " " << dRy << " " << dRz << " " << S_nonzero_number[count_n] << std::endl;
- }
+ g1[ispin] << dRx << " " << dRy << " " << dRz << " " << H_nonzero_num[ispin][count] << std::endl;
}
+ g2 << dRx << " " << dRy << " " << dRz << " " << S_nonzero_num[count] << std::endl;
+ }
+ }
- if (H_nonzero_number[count_n] == 0)
+ for (int ispin = 0; ispin < spin_loop; ++ispin)
+ {
+ if (H_nonzero_num[ispin][count] == 0)
+ {
+ // if (GlobalV::DRANK == 0)
+ // {
+ // if (!binary)
+ // {
+ // g1[ispin] << std::endl;
+ // g1[ispin] << std::endl;
+ // for (int index = 0; index < GlobalV::NLOCAL+1; ++index)
+ // {
+ // g1[ispin] << 0 << " ";
+ // }
+ // g1[ispin] << std::endl;
+ // }
+ // }
+ }
+ else
+ {
+ if (GlobalV::NSPIN != 4)
{
- // if (GlobalV::DRANK == 0)
- // {
- // if (!binary)
- // {
- // g1 << std::endl;
- // g1 << std::endl;
- // for (int index = 0; index < GlobalV::NLOCAL+1; ++index)
- // {
- // g1 << 0 << " ";
- // }
- // g1 << std::endl;
- // }
- // }
+ output_single_R(g1[ispin], HR_sparse_ptr[ispin][R_coor], sparse_threshold, binary);
}
else
{
- if (GlobalV::NSPIN != 4)
- {
- output_single_R(g1, HR_sparse_ptr[ix][iy][iz], sparse_threshold, binary);
- }
- else
- {
- output_soc_single_R(g1, HR_soc_sparse_ptr[ix][iy][iz], sparse_threshold, binary);
- }
+ output_soc_single_R(g1[ispin], HR_soc_sparse_ptr[R_coor], sparse_threshold, binary);
}
+ }
+ }
- if (S_nonzero_number[count_n] == 0)
+ if (S_nonzero_num[count] == 0)
+ {
+ // if (!binary)
+ // {
+ // if (GlobalV::DRANK == 0)
+ // {
+ // g2 << std::endl;
+ // g2 << std::endl;
+ // for (int index = 0; index < GlobalV::NLOCAL+1; ++index)
+ // {
+ // g2 << 0 << " ";
+ // }
+ // g2 << std::endl;
+ // }
+ // }
+ }
+ else
+ {
+ if (GlobalV::NSPIN != 4)
+ {
+ output_single_R(g2, SR_sparse_ptr[R_coor], sparse_threshold, binary);
+ }
+ else
+ {
+ output_soc_single_R(g2, SR_soc_sparse_ptr[R_coor], sparse_threshold, binary);
+ }
+ }
+
+ count++;
+
+ }
+
+ if(GlobalV::DRANK==0)
+ {
+ for (int ispin = 0; ispin < spin_loop; ++ispin) g1[ispin].close();
+ g2.close();
+ }
+
+ for (int ispin = 0; ispin < spin_loop; ++ispin)
+ {
+ delete[] H_nonzero_num[ispin];
+ H_nonzero_num[ispin] = nullptr;
+ }
+ delete[] S_nonzero_num;
+ S_nonzero_num = nullptr;
+
+ ModuleBase::timer::tick("HS_Matrix","save_HSR_sparse");
+ return;
+}
+
+void HS_Matrix::save_SR_sparse(
+ const double &sparse_threshold,
+ const bool &binary,
+ const std::string &SR_filename
+)
+{
+ ModuleBase::TITLE("HS_Matrix","save_SR_sparse");
+ ModuleBase::timer::tick("HS_Matrix","save_SR_sparse");
+
+ auto &all_R_coor_ptr = GlobalC::LM.all_R_coor;
+ auto &SR_sparse_ptr = GlobalC::LM.SR_sparse;
+ auto &SR_soc_sparse_ptr = GlobalC::LM.SR_soc_sparse;
+
+ int total_R_num = all_R_coor_ptr.size();
+ int output_R_number = 0;
+ int *S_nonzero_num = nullptr;
+
+ S_nonzero_num = new int[total_R_num];
+ ModuleBase::GlobalFunc::ZEROS(S_nonzero_num, total_R_num);
+
+ int count = 0;
+ for (auto &R_coor : all_R_coor_ptr)
+ {
+ if (GlobalV::NSPIN != 4)
+ {
+ auto iter = SR_sparse_ptr.find(R_coor);
+ if (iter != SR_sparse_ptr.end())
+ {
+ for (auto &row_loop : iter->second)
{
- // if (!binary)
- // {
- // if (GlobalV::DRANK == 0)
- // {
- // g2 << std::endl;
- // g2 << std::endl;
- // for (int index = 0; index < GlobalV::NLOCAL+1; ++index)
- // {
- // g2 << 0 << " ";
- // }
- // g2 << std::endl;
- // }
- // }
+ S_nonzero_num[count] += row_loop.second.size();
}
- else
+ }
+ }
+ else
+ {
+ auto iter = SR_soc_sparse_ptr.find(R_coor);
+ if (iter != SR_soc_sparse_ptr.end())
+ {
+ for (auto &row_loop : iter->second)
{
- if (GlobalV::NSPIN != 4)
- {
- output_single_R(g2, SR_sparse_ptr[ix][iy][iz], sparse_threshold, binary);
- }
- else
- {
- output_soc_single_R(g2, SR_soc_sparse_ptr[ix][iy][iz], sparse_threshold, binary);
- }
+ S_nonzero_num[count] += row_loop.second.size();
}
+ }
+ }
- count_n++;
+ count++;
+ }
+ Parallel_Reduce::reduce_int_all(S_nonzero_num, total_R_num);
+
+ for (int index = 0; index < total_R_num; ++index)
+ {
+ if (S_nonzero_num[index] != 0)
+ {
+ output_R_number++;
+ }
+ }
+
+ std::stringstream sss;
+ sss << SR_filename;
+ std::ofstream g2;
+
+ if(GlobalV::DRANK==0)
+ {
+ if (binary)
+ {
+ g2.open(sss.str().c_str(), ios::binary);
+ g2.write(reinterpret_cast(&GlobalV::NLOCAL), sizeof(int));
+ g2.write(reinterpret_cast(&output_R_number), sizeof(int));
+ }
+ else
+ {
+ g2.open(sss.str().c_str());
+ g2 << "Matrix Dimension of S(R): " << GlobalV::NLOCAL <(&dRx), sizeof(int));
+ g2.write(reinterpret_cast(&dRy), sizeof(int));
+ g2.write(reinterpret_cast(&dRz), sizeof(int));
+ g2.write(reinterpret_cast(&S_nonzero_num[count]), sizeof(int));
+ }
+ else
+ {
+ g2 << dRx << " " << dRy << " " << dRz << " " << S_nonzero_num[count] << std::endl;
}
}
+
+ if (GlobalV::NSPIN != 4)
+ {
+ output_single_R(g2, SR_sparse_ptr[R_coor], sparse_threshold, binary);
+ }
+ else
+ {
+ output_soc_single_R(g2, SR_soc_sparse_ptr[R_coor], sparse_threshold, binary);
+ }
+
+ count++;
+
}
if(GlobalV::DRANK==0)
{
- g1.close();
g2.close();
}
- delete[] H_nonzero_number;
- delete[] S_nonzero_number;
- H_nonzero_number = nullptr;
- S_nonzero_number = nullptr;
+ delete[] S_nonzero_num;
+ S_nonzero_num = nullptr;
- ModuleBase::timer::tick("HS_Matrix","save_HSR_sparse");
+ ModuleBase::timer::tick("HS_Matrix","save_SR_sparse");
return;
}
@@ -1105,9 +1301,10 @@ void HS_Matrix::output_single_R(std::ofstream &ofs, const std::map= 0)
@@ -1129,7 +1326,7 @@ void HS_Matrix::output_single_R(std::ofstream &ofs, const std::map sparse_threshold)
+ if (std::abs(line[col]) > sparse_threshold)
{
if (binary)
{
@@ -1151,11 +1348,14 @@ void HS_Matrix::output_single_R(std::ofstream &ofs, const std::map[GlobalV::NLOCAL];
for(int row = 0; row < GlobalV::NLOCAL; ++row)
{
- line = new std::complex[GlobalV::NLOCAL];
+ // line = new std::complex[GlobalV::NLOCAL];
ModuleBase::GlobalFunc::ZEROS(line, GlobalV::NLOCAL);
if(GlobalC::ParaO.trace_loc_row[row] >= 0)
@@ -1238,7 +1439,7 @@ void HS_Matrix::output_soc_single_R(std::ofstream &ofs, const std::map sparse_threshold)
+ if (std::abs(line[col]) > sparse_threshold)
{
if (binary)
{
@@ -1261,11 +1462,14 @@ void HS_Matrix::output_soc_single_R(std::ofstream &ofs, const std::map *H, const std::complex *S, bool bit);
-
- void save_HSR_tr(const int current_spin); //LiuXh add 2019-07-15
-
- // jingan add 2021-6-4
- void save_HSR_sparse(const int ¤t_spin, const double &sparse_threshold, const bool &binary);
- void output_single_R(std::ofstream &ofs, const std::map> &XR, const double &sparse_threshold, const bool &binary);
- void output_soc_single_R(std::ofstream &ofs, const std::map>> &XR, const double &sparse_threshold, const bool &binary);
+ void saving_HS(const double *Hloc, const double* Sloc, bool bit, const int &out_hs);
+
+ void save_HS(const double *H, const double *S, bool bit);
+
+ void save_HS_complex(const std::complex *H, const std::complex *S, bool bit);
+
+ void save_HSR_tr(const int current_spin); //LiuXh add 2019-07-15
+
+ // jingan add 2021-6-4, modify 2021-12-2
+ void save_HSR_sparse(
+ const double &sparse_threshold,
+ const bool &binary,
+ const std::string &SR_filename,
+ const std::string &HR_filename_up,
+ const std::string &HR_filename_down
+ );
+ void save_SR_sparse(
+ const double &sparse_threshold,
+ const bool &binary,
+ const std::string &SR_filename
+ );
+ void output_single_R(std::ofstream &ofs, const std::map> &XR, const double &sparse_threshold, const bool &binary);
+ void output_soc_single_R(std::ofstream &ofs, const std::map>> &XR, const double &sparse_threshold, const bool &binary);
// mohan comment out 2021-02-10
-// void save_HS_ccf(const int &iter, const int &Hnnz, const int *colptr_H, const int *rowind_H,
-// const double *nzval_H, const double *nzval_S, bool bit);
+// void save_HS_ccf(const int &iter, const int &Hnnz, const int *colptr_H, const int *rowind_H,
+// const double *nzval_H, const double *nzval_S, bool bit);
- void saving_HS_complex(std::complex *Hloc, std::complex* Sloc, bool bit, const int &out_hs); //LiuXh, 2017-03-21
+ void saving_HS_complex(std::complex *Hloc, std::complex* Sloc, bool bit, const int &out_hs); //LiuXh, 2017-03-21
- void save_HS_complex(std::complex *H, std::complex *S, bool bit); //LiuXh, 2017-03-21
+ void save_HS_complex(std::complex *H, std::complex *S, bool bit); //LiuXh, 2017-03-21
}
#endif
diff --git a/source/src_io/write_HS_R.cpp b/source/src_io/write_HS_R.cpp
index b4edf71e0fa..086c2165027 100644
--- a/source/src_io/write_HS_R.cpp
+++ b/source/src_io/write_HS_R.cpp
@@ -3,23 +3,27 @@
#include "../src_pw/global.h"
#include "write_HS.h"
-
-void LOOP_ions::output_HS_R(void)
+// if 'binary=true', output binary file.
+// The 'sparse_threshold' is the accuracy of the sparse matrix.
+// If the absolute value of the matrix element is less than or equal to the 'sparse_threshold', it will be ignored.
+void LOOP_ions::output_HS_R(
+ const std::string &SR_filename,
+ const std::string &HR_filename_up,
+ const std::string HR_filename_down,
+ const bool &binary,
+ const double &sparse_threshold
+)
{
ModuleBase::TITLE("LOOP_ions","output_HS_R");
ModuleBase::timer::tick("LOOP_ions","output_HS_R");
-
- // add by jingan for out r_R matrix 2019.8.14
- if(INPUT.out_r_matrix)
- {
- cal_r_overlap_R r_matrix;
- r_matrix.init();
- r_matrix.out_r_overlap_R(GlobalV::NSPIN);
- }
-
- // Parameters for HR and SR output
- double sparse_threshold = 1e-10;
- bool binary = false; // output binary file
+
+ // add by jingan for out r_R matrix 2019.8.14
+ if(INPUT.out_r_matrix)
+ {
+ cal_r_overlap_R r_matrix;
+ r_matrix.init();
+ r_matrix.out_r_overlap_R(GlobalV::NSPIN);
+ }
if(GlobalV::NSPIN==1||GlobalV::NSPIN==4)
{
@@ -28,10 +32,8 @@ void LOOP_ions::output_HS_R(void)
// GlobalC::UHM.GK.distribute_pvpR_tr();
// HS_Matrix::save_HSR_tr(0);
- // jingan add 2021-6-4
+ // jingan add 2021-6-4, modify 2021-12-2
GlobalC::UHM.calculate_HSR_sparse(0, sparse_threshold);
- HS_Matrix::save_HSR_sparse(0, sparse_threshold, binary);
- GlobalC::UHM.destroy_all_HSR_sparse();
}
///*
else if(GlobalV::NSPIN==2)
@@ -46,13 +48,13 @@ void LOOP_ions::output_HS_R(void)
// {
// GlobalC::pot.vr_eff1[ir] = GlobalC::pot.vr_eff( GlobalV::CURRENT_SPIN, ir);
// }
-
+
// if(!GlobalV::GAMMA_ONLY_LOCAL)
// {
// if(GlobalV::VL_IN_H)
// {
- // //GlobalC::UHM.GK.cal_vlocal_k(GlobalC::pot.vrs1,GridT);
- // GlobalC::UHM.GK.cal_vlocal_k(GlobalC::pot.vr_eff1, GlobalC::GridT, GlobalV::CURRENT_SPIN);
+ // //GlobalC::UHM.GK.cal_vlocal_k(GlobalC::pot.vrs1,GridT);
+ // GlobalC::UHM.GK.cal_vlocal_k(GlobalC::pot.vr_eff1, GlobalC::GridT, GlobalV::CURRENT_SPIN);
// }
// }
// GlobalC::UHM.GK.cal_vlocal_R(GlobalV::CURRENT_SPIN);
@@ -75,22 +77,24 @@ void LOOP_ions::output_HS_R(void)
{
GlobalC::pot.vr_eff1[ir] = GlobalC::pot.vr_eff( GlobalV::CURRENT_SPIN, ir);
}
-
+
if(!GlobalV::GAMMA_ONLY_LOCAL)
{
if(GlobalV::VL_IN_H)
{
- //GlobalC::UHM.GK.cal_vlocal_k(GlobalC::pot.vrs1,GridT);
- GlobalC::UHM.GK.cal_vlocal_k(GlobalC::pot.vr_eff1, GlobalC::GridT, GlobalV::CURRENT_SPIN);
+ //GlobalC::UHM.GK.cal_vlocal_k(GlobalC::pot.vrs1,GridT);
+ GlobalC::UHM.GK.cal_vlocal_k(GlobalC::pot.vr_eff1, GlobalC::GridT, GlobalV::CURRENT_SPIN);
}
}
+
GlobalC::UHM.calculate_HSR_sparse(GlobalV::CURRENT_SPIN, sparse_threshold);
- HS_Matrix::save_HSR_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, binary);
- GlobalC::UHM.destroy_all_HSR_sparse();
}
}
}
+ HS_Matrix::save_HSR_sparse(sparse_threshold, binary, SR_filename, HR_filename_up, HR_filename_down);
+ GlobalC::UHM.destroy_all_HSR_sparse();
+
if(!GlobalV::GAMMA_ONLY_LOCAL) //LiuXh 20181011
{
GlobalC::UHM.GK.destroy_pvpR();
@@ -99,3 +103,17 @@ void LOOP_ions::output_HS_R(void)
ModuleBase::timer::tick("LOOP_ions","output_HS_R");
return;
}
+
+
+void LOOP_ions::output_SR(const std::string &SR_filename, const bool &binary, const double &sparse_threshold)
+{
+ ModuleBase::TITLE("LOOP_ions","output_SR");
+ ModuleBase::timer::tick("LOOP_ions","output_SR");
+
+ GlobalC::UHM.calculate_SR_sparse(sparse_threshold);
+ HS_Matrix::save_SR_sparse(sparse_threshold, binary, SR_filename);
+ GlobalC::UHM.destroy_all_HSR_sparse();
+
+ ModuleBase::timer::tick("LOOP_ions","output_SR");
+ return;
+}
\ No newline at end of file
diff --git a/source/src_io/write_input.cpp b/source/src_io/write_input.cpp
index 74a548c132f..f824400049e 100644
--- a/source/src_io/write_input.cpp
+++ b/source/src_io/write_input.cpp
@@ -59,6 +59,7 @@ void Input::Print(const std::string &fn)const
ModuleBase::GlobalFunc::OUTP(ofs,"out_charge",out_charge,">0 output charge density for selected electron steps");
ModuleBase::GlobalFunc::OUTP(ofs,"out_potential",out_potential,"output realspace potential");
ModuleBase::GlobalFunc::OUTP(ofs,"out_wf",out_wf,"output wave functions");
+ ModuleBase::GlobalFunc::OUTP(ofs,"out_wf_r",out_wf_r,"output wave functions in realspace");
ModuleBase::GlobalFunc::OUTP(ofs,"out_dos",out_dos,"output energy and dos");
ModuleBase::GlobalFunc::OUTP(ofs,"out_band",out_band,"output energy and band structure");
ModuleBase::GlobalFunc::OUTP(ofs,"restart_save",restart_save,"print to disk every step for restart");
@@ -174,6 +175,13 @@ void Input::Print(const std::string &fn)const
ModuleBase::GlobalFunc::OUTP(ofs,"rcut_lj",mdp.rcut_lj,"cutoff radius of LJ potential");
ModuleBase::GlobalFunc::OUTP(ofs,"epsilon_lj",mdp.epsilon_lj,"the value of epsilon for LJ potential");
ModuleBase::GlobalFunc::OUTP(ofs,"sigma_lj",mdp.sigma_lj,"the value of sigma for LJ potential");
+ ModuleBase::GlobalFunc::OUTP(ofs,"direction",mdp.direction,"the direction of shock wave");
+ ModuleBase::GlobalFunc::OUTP(ofs,"velocity",mdp.velocity,"the velocity of shock wave");
+ ModuleBase::GlobalFunc::OUTP(ofs,"viscosity",mdp.viscosity,"artificial viscosity");
+ ModuleBase::GlobalFunc::OUTP(ofs,"tscale",mdp.tscale,"reduction in initial temperature");
+ ModuleBase::GlobalFunc::OUTP(ofs,"md_tfreq",mdp.tfreq,"oscillation frequency, used to determine Qmass of NHC");
+ ModuleBase::GlobalFunc::OUTP(ofs,"md_damp",mdp.damp,"damping parameter (time units) used to add force in Langevin method");
+
ofs << "\n#Parameters (11.Efield)" << std::endl;
ModuleBase::GlobalFunc::OUTP(ofs,"efield",efield,"add electric field");
@@ -197,7 +205,6 @@ void Input::Print(const std::string &fn)const
ModuleBase::GlobalFunc::OUTP(ofs,"test_stress", test_stress, "test the force");
ofs << "\n#Parameters (13.Other Methods)" << std::endl;
- ModuleBase::GlobalFunc::OUTP(ofs,"mlwf_flag",mlwf_flag,"turn MLWF on or off");
ModuleBase::GlobalFunc::OUTP(ofs,"opt_epsilon2",opt_epsilon2,"calculate the dielectic function");
ModuleBase::GlobalFunc::OUTP(ofs,"opt_nbands",opt_nbands,"number of bands for optical calculation");
diff --git a/source/src_io/write_rho_cube.cpp b/source/src_io/write_rho_cube.cpp
index 8d491e8bf77..ff7ebc375cf 100644
--- a/source/src_io/write_rho_cube.cpp
+++ b/source/src_io/write_rho_cube.cpp
@@ -110,6 +110,7 @@ void Charge::write_rho_cube(
{
start_z[ip] = start_z[ip-1]+num_z[ip-1];
}
+ delete[] num_z;
// which_ip: found iz belongs to which ip.
int *which_ip = new int[GlobalC::pw.ncz];
diff --git a/source/src_io/write_wfc_realspace.cpp b/source/src_io/write_wfc_realspace.cpp
new file mode 100644
index 00000000000..2f8fe4b14b3
--- /dev/null
+++ b/source/src_io/write_wfc_realspace.cpp
@@ -0,0 +1,160 @@
+//======================
+// AUTHOR : Peize Lin
+// DATE : 2021-11-21
+//======================
+
+#include "write_wfc_realspace.h"
+#include "src_pw/global.h"
+#include "module_base/tool_title.h"
+#include
+#include
+#include
+
+namespace Write_Wfc_Realspace
+{
+ // write ||wfc_r|| for all k-points and all bands
+ // Input: wfc_g[ik](ib,ig)
+ // loop order is for(z){for(y){for(x)}}
+ void write_wfc_realspace_1(const ModuleBase::ComplexMatrix*const wfc_g, const std::string &folder_name)
+ {
+ ModuleBase::TITLE("Write_Wfc_Realspace", "write_wfc_realspace_1");
+
+ const string outdir = GlobalV::global_out_dir + folder_name + "/";
+ const std::string command0 = "test -d " + outdir + " || mkdir " + outdir;
+ if(GlobalV::MY_RANK==0)
+ system( command0.c_str() );
+
+#ifdef __MPI
+ std::vector mpi_requests;
+#endif
+ for(int ik=0; ik> wfc_r = cal_wfc_r(wfc_g[ik], ik, ib);
+
+ std::vector wfc_r2(wfc_r.size());
+ for(int ir=0; ir
+ // rank0 k0 k1 k2 k3 k4 k5
+ // \ \ \ \ \ \
+ // rank1 k0 k1 k2 k3 k4 k5
+ // \ \ \ \ \ \
+ // rank2 k0 k1 k2 k3 k4 k5
+
+
+
+ // Input: wfc_g(ib,ig)
+ // Output: wfc_r[ir]
+ std::vector> cal_wfc_r(const ModuleBase::ComplexMatrix &wfc_g, const int ik, const int ib)
+ {
+ ModuleBase::GlobalFunc::ZEROS(GlobalC::UFFT.porter, GlobalC::pw.nrxx);
+ std::vector> wfc_r(GlobalC::pw.nrxx);
+ for(int ig=0; ig &chg_r, const std::string &file_name, MPI_Request &mpi_request)
+#else
+ void write_charge_realspace_1(const std::vector &chg_r, const std::string &file_name)
+#endif
+ {
+ std::ofstream ofs;
+
+#ifdef __MPI
+ constexpr int mpi_tag=100;
+ if(GlobalV::RANK_IN_POOL==0)
+ {
+#endif
+ ofs.open(file_name);
+
+ ofs<<"calculated by ABACUS"<
+#include
+#include
+
+#ifdef __MPI
+#include
+#endif
+
+namespace Write_Wfc_Realspace
+{
+ // write ||wfc_r|| for all k-points and all bands
+ // Input: wfc_g[ik](ib,ig)
+ // loop order is for(z){for(y){for(x)}}
+ void write_wfc_realspace_1(const ModuleBase::ComplexMatrix*const wfc_g, const std::string &folder_name);
+
+ // Input: wfc_g(ib,ig)
+ // Output: wfc_r[ir]
+ std::vector