diff --git a/README.md b/README.md index 3295bff03..5d8214b27 100644 --- a/README.md +++ b/README.md @@ -82,6 +82,7 @@ The `System` or `LabeledSystem` can be constructed from the following file forma | Gromacs | gro | True | False | System | 'gromacs/gro' | | ABACUS | STRU | False | True | LabeledSystem | 'abacus/scf' | | ABACUS | cif | True | True | LabeledSystem | 'abacus/md' | +| ABACUS | STRU | True | True | LabeledSystem | 'abacus/relax' | | ase | structure | True | True | MultiSystems | 'ase/structure' | diff --git a/dpdata/abacus/relax.py b/dpdata/abacus/relax.py new file mode 100644 index 000000000..3b334ab92 --- /dev/null +++ b/dpdata/abacus/relax.py @@ -0,0 +1,139 @@ +import os,sys +import numpy as np +from .scf import bohr2ang, kbar2evperang3, get_geometry_in, get_cell, get_coords + +# Read in geometries from an ABACUS RELAX(CELL-RELAX) trajectory in OUT.XXXX/runnning_relax/cell-relax.log. + +def get_log_file(fname, inlines): + suffix = "ABACUS" + calculation = "scf" + for line in inlines: + if "suffix" in line and "suffix"==line.split()[0]: + suffix = line.split()[1] + elif "calculation" in line and "calculation" == line.split()[0]: + calculation = line.split()[1] + logf = os.path.join(fname, "OUT.%s/running_%s.log"%(suffix,calculation)) + return logf + +def get_coords_from_log(loglines,natoms): + ''' + NOTICE: unit of coords and cells is Angstrom + ''' + natoms_log = 0 + for line in loglines: + if line[13:41] == "number of atom for this type": + natoms_log += int(line.split()[-1]) + + assert(natoms_log>0 and natoms_log == natoms),"ERROR: detected atom number in log file is %d" % natoms + + energy = [] + cells = [] + coords = [] + force = [] + stress = [] + + for i in range(len(loglines)): + line = loglines[i] + if line[18:41] == "lattice constant (Bohr)": + a0 = float(line.split()[-1]) + elif len(loglines[i].split()) >=2 and loglines[i].split()[1] == 'COORDINATES': + coords.append([]) + direct_coord = False + if loglines[i].split()[0] == 'DIRECT': + direct_coord = True + for k in range(2,2+natoms): + coords[-1].append(list(map(lambda x: float(x),loglines[i+k].split()[1:4]))) + elif loglines[i].split()[0] == 'CARTESIAN': + for k in range(2,2+natoms): + coords[-1].append(list(map(lambda x: float(x)*a0,loglines[i+k].split()[1:4]))) + else: + assert(False),"Unrecongnized coordinate type, %s, line:%d" % (loglines[i].split()[0],i) + + converg = True + for j in range(i): + if loglines[i-j-1][1:36] == 'Ion relaxation is not converged yet': + converg = False + break + elif loglines[i-j-1][1:29] == 'Ion relaxation is converged!': + converg = True + break + + if converg: + for j in range(i+1,len(loglines)): + if loglines[j][1:56] == "Lattice vectors: (Cartesian coordinate: in unit of a_0)": + cells.append([]) + for k in range(1,4): + cells[-1].append(list(map(lambda x:float(x)*a0,loglines[j+k].split()[0:3]))) + break + else: + cells.append(cells[-1]) + + if direct_coord: + coords[-1] = coords[-1].dot(cells[-1]) + + elif line[4:15] == "TOTAL-FORCE": + force.append([]) + for j in range(5,5+natoms): + force[-1].append(list(map(lambda x:float(x),loglines[i+j].split()[1:4]))) + elif line[1:13] == "TOTAL-STRESS": + stress.append([]) + for j in range(4,7): + stress[-1].append(list(map(lambda x:float(x),loglines[i+j].split()[0:3]))) + elif line[1:14] == "final etot is": + energy.append(float(line.split()[-2])) + + assert(len(cells) == len(coords) or len(cells)+1 == len(coords)),"ERROR: detected %d coordinates and %d cells" % (len(coords),len(cells)) + if len(cells)+1 == len(coords): del(coords[-1]) + + energy = np.array(energy) + cells = np.array(cells) + coords = np.array(coords) + stress = np.array(stress) + force = np.array(force) + + cells *= bohr2ang + coords *= bohr2ang + + virial = np.zeros([len(cells), 3, 3]) + for i in range(len(cells)): + volume = np.linalg.det(cells[i, :, :].reshape([3, 3])) + virial[i] = stress[i] * kbar2evperang3 * volume + + return energy,cells,coords,force,stress,virial + +def get_frame (fname): + if type(fname) == str: + # if the input parameter is only one string, it is assumed that it is the + # base directory containing INPUT file; + path_in = os.path.join(fname, "INPUT") + else: + raise RuntimeError('invalid input') + with open(path_in, 'r') as fp: + inlines = fp.read().split('\n') + geometry_path_in = get_geometry_in(fname, inlines) # base dir of STRU + with open(geometry_path_in, 'r') as fp: + geometry_inlines = fp.read().split('\n') + celldm, cell = get_cell(geometry_inlines) + atom_names, natoms, types, coord_tmp = get_coords(celldm, cell, geometry_inlines, inlines) + + logf = get_log_file(fname, inlines) + assert(os.path.isfile(logf)),"Error: can not find %s" % logf + with open(logf) as f1: lines = f1.readlines() + + atomnumber = 0 + for i in natoms: atomnumber += i + energy,cells,coords,force,stress,virial = get_coords_from_log(lines,atomnumber) + + data = {} + data['atom_names'] = atom_names + data['atom_numbs'] = natoms + data['atom_types'] = types + data['cells'] = cells + data['coords'] = coords + data['energies'] = energy + data['forces'] = force + data['virials'] = virial + data['stress'] = stress + data['orig'] = np.zeros(3) + + return data diff --git a/dpdata/plugins/abacus.py b/dpdata/plugins/abacus.py index e4e1c291e..c219053b7 100644 --- a/dpdata/plugins/abacus.py +++ b/dpdata/plugins/abacus.py @@ -1,5 +1,6 @@ import dpdata.abacus.scf import dpdata.abacus.md +import dpdata.abacus.relax from dpdata.format import Format @Format.register("abacus/stru") @@ -51,3 +52,11 @@ class AbacusMDFormat(Format): #@Format.post("rot_lower_triangular") def from_labeled_system(self, file_name, **kwargs): return dpdata.abacus.md.get_frame(file_name) + +@Format.register("abacus/relax") +@Format.register("abacus/pw/relax") +@Format.register("abacus/lcao/relax") +class AbacusRelaxFormat(Format): + #@Format.post("rot_lower_triangular") + def from_labeled_system(self, file_name, **kwargs): + return dpdata.abacus.relax.get_frame(file_name) diff --git a/dpdata/system.py b/dpdata/system.py index a919d31fc..3c9ccec8c 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -195,7 +195,8 @@ def __init__ (self, - ``qe/cp/traj``: Quantum Espresso CP trajectory files. should have: file_name+'.in' and file_name+'.pos' - ``qe/pw/scf``: Quantum Espresso PW single point calculations. Both input and output files are required. If file_name is a string, it denotes the output file name. Input file name is obtained by replacing 'out' by 'in' from file_name. Or file_name is a list, with the first element being the input file name and the second element being the output filename. - ``abacus/scf``: ABACUS pw/lcao scf. The directory containing INPUT file is required. - - ``abacus/md``: ABACUS pw/lcao MD. The directory containing INPUT file is required. + - ``abacus/md``: ABACUS pw/lcao MD. The directory containing INPUT file is required. + - ``abacus/relax``: ABACUS pw/lcao relax or cell-relax. The directory containing INPUT file is required. - ``siesta/output``: siesta SCF output file - ``siesta/aimd_output``: siesta aimd output file - ``pwmat/atom.config``: pwmat atom.config diff --git a/tests/abacus.relax/INPUT b/tests/abacus.relax/INPUT new file mode 100644 index 000000000..46ab55eac --- /dev/null +++ b/tests/abacus.relax/INPUT @@ -0,0 +1,29 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix abacus +calculation cell-relax +ntype 2 +nbands 6 +symmetry 1 + +#Parameters (2.Iteration) +ecutwfc 50 +scf_thr 1e-7 +scf_nmax 50 + +relax_nmax 5 + +#Parameters (3.Basis) +basis_type pw + +#Parameters (4.Smearing) +smearing_method gaussian +smearing_sigma 0.02 + +#Parameters (5.Mixing) +mixing_type pulay +mixing_beta 0.4 + +#Parameters (6.Deepks) +cal_force 1 +cal_stress 1 diff --git a/tests/abacus.relax/KPT b/tests/abacus.relax/KPT new file mode 100644 index 000000000..c289c0158 --- /dev/null +++ b/tests/abacus.relax/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/abacus.relax/OUT.abacus/INPUT b/tests/abacus.relax/OUT.abacus/INPUT new file mode 100644 index 000000000..88e098762 --- /dev/null +++ b/tests/abacus.relax/OUT.abacus/INPUT @@ -0,0 +1,250 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix abacus #the name of main output directory +latname test #the name of lattice name +stru_file STRU #the filename of file containing atom positions +kpoint_file KPT #the name of file containing k points +pseudo_dir #the directory containing pseudo files +orbital_dir #the directory containing orbital files +pseudo_type auto #the type pseudo files +pseudo_rcut 15 #cut-off radius for radial integration +pseudo_mesh 0 #0: use our own mesh to do radial renormalization; 1: use mesh as in QE +lmaxmax 2 #maximum of l channels used +dft_functional default #exchange correlation functional +calculation cell-relax #test; scf; relax; nscf; ienvelope; istate; sto-scf; sto-md +ntype 2 #atom species number +nspin 1 #1: single spin; 2: up and down spin; 4: noncollinear spin +kspacing 0 #unit in 1/bohr, should be > 0, default is 0 which means read KPT file +nbands 6 #number of bands +nbands_sto 256 #number of stochastic bands +nbands_istate 5 #number of bands around Fermi level for istate calulation +nche_sto 100 #number of orders for Chebyshev expansion in stochastic DFT +symmetry 1 #turn symmetry on or off +init_vel 0 #read velocity from STRU or not +symmetry_prec 1e-05 #accuracy for symmetry +nelec 0 #input number of electrons +tot_magnetization 0 #total magnetization of the system +out_mul 0 # mulliken charge or not +noncolin 0 #using non-collinear-spin +lspinorb 0 #consider the spin-orbit interaction +kpar 1 #devide all processors into kpar groups and k points will be distributed among each group +bndpar 1 #devide all processors into bndpar groups and bands will be distributed among each group + +#Parameters (2.PW) +ecutwfc 50 ##energy cutoff for wave functions +pw_diag_nmax 50 #max iteration number for cg +diago_cg_prec 1 #diago_cg_prec +pw_diag_thr 0.01 #threshold for eigenvalues is cg electron iterations +scf_thr 1e-07 #charge density error +init_wfc atomic #start wave functions are from 'atomic', 'atomic+random', 'random' or 'file' +init_chg atomic #start charge is from 'atomic' or file +chg_extrap atomic #atomic; first-order; second-order; dm:coefficients of SIA +out_chg 0 #>0 output charge density for selected electron steps +out_pot 0 #output realspace potential +out_wfc_pw 0 #output wave functions +out_wfc_r 0 #output wave functions in realspace +out_dos 0 #output energy and dos +out_band 0 #output energy and band structure +out_proj_band 0 #output projected band structure +restart_save 0 #print to disk every step for restart +restart_load 0 #restart from disk +read_file_dir auto #directory of files for reading +nx 0 #number of points along x axis for FFT grid +ny 0 #number of points along y axis for FFT grid +nz 0 #number of points along z axis for FFT grid +cell_factor 1.2 #used in the construction of the pseudopotential tables + +#Parameters (3.Stochastic DFT) +method_sto 1 #1: slow and save memory, 2: fast and waste memory +nbands_sto 256 #number of stochstic orbitals +nche_sto 100 #Chebyshev expansion orders +emin_sto 0 #trial energy to guess the lower bound of eigen energies of the Hamitonian operator +emax_sto 0 #trial energy to guess the upper bound of eigen energies of the Hamitonian operator +seed_sto 0 #the random seed to generate stochastic orbitals +initsto_freq 1000 #frequency to generate new stochastic orbitals when running md + +#Parameters (4.Relaxation) +ks_solver cg #cg; dav; lapack; genelpa; hpseps; scalapack_gvx; cusolver +scf_nmax 50 ##number of electron iterations +out_force 0 #output the out_force or not +relax_nmax 5 #number of ion iteration steps +out_stru 0 #output the structure files after each ion step +force_thr 0.001 #force threshold, unit: Ry/Bohr +force_thr_ev 0.0257112 #force threshold, unit: eV/Angstrom +force_thr_ev2 0 #force invalid threshold, unit: eV/Angstrom +relax_cg_thr 0.5 #threshold for switching from cg to bfgs, unit: eV/Angstrom +stress_thr 0.01 #stress threshold +press1 0 #target pressure, unit: KBar +press2 0 #target pressure, unit: KBar +press3 0 #target pressure, unit: KBar +relax_bfgs_w1 0.01 #wolfe condition 1 for bfgs +relax_bfgs_w2 0.5 #wolfe condition 2 for bfgs +relax_bfgs_rmax 0.8 #maximal trust radius, unit: Bohr +relax_bfgs_rmin 1e-05 #minimal trust radius, unit: Bohr +relax_bfgs_init 0.5 #initial trust radius, unit: Bohr +cal_stress 1 #calculate the stress or not +fixed_axes None #which axes are fixed +relax_method cg #bfgs; sd; cg; cg_bfgs; +out_level ie #ie(for electrons); i(for ions); +out_dm 0 #>0 output density matrix +deepks_out_labels 0 #>0 compute descriptor for deepks +deepks_scf 0 #>0 add V_delta to Hamiltonian +deepks_bandgap 0 #>0 for bandgap label +deepks_out_unittest 0 #if set 1, prints intermediate quantities that shall be used for making unit test +deepks_model #file dir of traced pytorch model: 'model.ptg +deepks_descriptor_lmax2 #lmax used in generating descriptor +deepks_descriptor_rcut0 #rcut used in generating descriptor +deepks_descriptor_ecut0 #ecut used in generating descriptor + +#Parameters (5.LCAO) +basis_type pw #PW; LCAO in pw; LCAO +search_radius -1 #input search radius (Bohr) +search_pbc 1 #input periodic boundary condition +lcao_ecut 0 #energy cutoff for LCAO +lcao_dk 0.01 #delta k for 1D integration in LCAO +lcao_dr 0.01 #delta r for 1D integration in LCAO +lcao_rmax 30 #max R for 1D two-center integration table +out_mat_hs 0 #output H and S matrix +out_mat_hs2 0 #output H(R) and S(R) matrix +out_mat_r 0 #output r(R) matrix +out_wfc_lcao 0 #ouput LCAO wave functions +bx 1 #division of an element grid in FFT grid along x +by 1 #division of an element grid in FFT grid along y +bz 1 #division of an element grid in FFT grid along z + +#Parameters (6.Smearing) +smearing_method gaussian #type of smearing_method: gauss; fd; fixed; mp; mp2; mv +smearing_sigma 0.02 #energy range for smearing + +#Parameters (7.Charge Mixing) +mixing_type pulay #plain; kerker; pulay; pulay-kerker; broyden +mixing_beta 0.4 #mixing parameter: 0 means no new charge +mixing_ndim 8 #mixing dimension in pulay +mixing_gg0 0 #mixing parameter in kerker + +#Parameters (8.DOS) +dos_emin_ev -15 #minimal range for dos +dos_emax_ev 15 #maximal range for dos +dos_edelta_ev 0.01 #delta energy for dos +dos_scale 0.01 #scale dos range by +dos_sigma 0.07 #gauss b coefficeinet(default=0.07) + +#Parameters (9.Molecular dynamics) +md_type 1 #choose ensemble +md_nstep 10 #md steps +md_ensolver FP #choose potential +md_dt 1 #time step +md_mnhc 4 #number of Nose-Hoover chains +md_tfirst -1 #temperature first +md_tlast -1 #temperature last +md_dumpfreq 1 #The period to dump MD information +md_restartfreq 5 #The period to output MD restart information +md_seed -1 #random seed for MD +md_restart 0 #whether restart +lj_rcut 8.5 #cutoff radius of LJ potential +lj_epsilon 0.01032 #the value of epsilon for LJ potential +lj_sigma 3.405 #the value of sigma for LJ potential +msst_direction 2 #the direction of shock wave +msst_vel 0 #the velocity of shock wave +msst_vis 0 #artificial viscosity +msst_tscale 0.01 #reduction in initial temperature +msst_qmass -1 #mass of thermostat +md_tfreq 0 #oscillation frequency, used to determine qmass of NHC +md_damp 1 #damping parameter (time units) used to add force in Langevin method + +#Parameters (10.Electric field and dipole correction) +efield_flag 0 #add electric field +dip_cor_flag 0 #dipole correction +efield_dir 2 #the direction of the electric field or dipole correction +efield_pos_max 0.5 #position of the maximum of the saw-like potential along crystal axis efield_dir +efield_pos_dec 0.1 #zone in the unit cell where the saw-like potential decreases +efield_amp 0 #amplitude of the electric field + +#Parameters (11.Test) +out_alllog 0 #output information for each processor, when parallel +nurse 0 #for coders +colour 0 #for coders, make their live colourful +t_in_h 1 #calculate the kinetic energy or not +vl_in_h 1 #calculate the local potential or not +vnl_in_h 1 #calculate the nonlocal potential or not +vh_in_h 1 #calculate the hartree potential or not +vion_in_h 1 #calculate the local ionic potential or not +test_force 0 #test the force +test_stress 0 #test the force + +#Parameters (13.vdW Correction) +vdw_method none #the method of calculating vdw (none ; d2 ; d3_0 ; d3_bj +vdw_s6 default #scale parameter of d2/d3_0/d3_bj +vdw_s8 default #scale parameter of d3_0/d3_bj +vdw_a1 default #damping parameter of d3_0/d3_bj +vdw_a2 default #damping parameter of d3_bj +vdw_d 20 #damping parameter of d2 +vdw_abc 0 #third-order term? +vdw_C6_file default #filename of C6 +vdw_C6_unit Jnm6/mol #unit of C6, Jnm6/mol or eVA6 +vdw_R0_file default #filename of R0 +vdw_R0_unit A #unit of R0, A or Bohr +vdw_model radius #expression model of periodic structure, radius or period +vdw_radius default #radius cutoff for periodic structure +vdw_radius_unit Bohr #unit of radius cutoff for periodic structure +vdw_cn_thr 40 #radius cutoff for cn +vdw_cn_thr_unit Bohr #unit of cn_thr, Bohr or Angstrom +vdw_period 3 3 3 #periods of periodic structure + +#Parameters (14.exx) +dft_functional default #no, hf, pbe0, hse or opt_orb +exx_hybrid_alpha 0.25 # +exx_hse_omega 0.11 # +exx_separate_loop 1 #0 or 1 +exx_hybrid_step 100 # +exx_lambda 0.3 # +exx_pca_threshold 0 # +exx_c_threshold 0 # +exx_v_threshold 0 # +exx_dm_threshold 0 # +exx_schwarz_threshold0 # +exx_cauchy_threshold0 # +exx_ccp_threshold 1e-08 # +exx_ccp_rmesh_times 10 # +exx_distribute_type htime #htime or kmeans1 or kmeans2 +exx_opt_orb_lmax 0 # +exx_opt_orb_ecut 0 # +exx_opt_orb_tolerence0 # + +#Parameters (16.tddft) +tddft 0 #calculate tddft or not +td_scf_thr 1e-09 #threshold for electronic iteration of tddft +td_dt 0.02 #time of ion step +td_force_dt 0.02 #time of force change +td_val_elec_01 1 #td_val_elec_01 +td_val_elec_02 1 #td_val_elec_02 +td_val_elec_03 1 #td_val_elec_03 +td_vext 0 #add extern potential or not +td_vext_dire 1 #extern potential direction +td_timescale 0.5 #extern potential td_timescale +td_vexttype 1 #extern potential type +td_vextout 0 #output extern potential or not +td_dipoleout 0 #output dipole or not +ocp 0 #change occupation or not +ocp_set none #set occupation + +#Parameters (17.berry_wannier) +berry_phase 0 #calculate berry phase or not +gdir 3 #calculate the polarization in the direction of the lattice vector +towannier90 0 #use wannier90 code interface or not +nnkpfile seedname.nnkp #the wannier90 code nnkp file name +wannier_spin up #calculate spin in wannier90 code interface + +#Parameters (18.implicit_solvation) +imp_sol 0 #calculate implicit solvation correction or not +eb_k 80 #the relative permittivity of the bulk solvent +tau 1.0798e-05 #the effective surface tension parameter +sigma_k 0.6 # the width of the diffuse cavity +nc_k 0.00037 # the cut-off charge density + +#Parameters (19.compensating_charge) +comp_chg 0 # add compensating charge +comp_q 0 # total charge of compensating charge +comp_l 1 # total length of compensating charge +comp_center 0 # center of compensating charge on dim +comp_dim 2 # dimension of compensating charge(x, y or z) diff --git a/tests/abacus.relax/OUT.abacus/STRU_ION_D b/tests/abacus.relax/OUT.abacus/STRU_ION_D new file mode 100644 index 000000000..80feeb742 --- /dev/null +++ b/tests/abacus.relax/OUT.abacus/STRU_ION_D @@ -0,0 +1,25 @@ +ATOMIC_SPECIES +H 1.008 ../potential/H_ONCV_PBE-1.0.upf +O 15.9994 ../potential/O_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1 + +LATTICE_VECTORS +28 0 0 #latvec1 +0 28 0 #latvec2 +0 0 28 #latvec3 + +ATOMIC_POSITIONS +Direct + +H #label +0 #magnetism +2 #number of atoms +0.569506832635 0.670672815952 0.298608009885 m 1 1 1 +0.492739894831 0.735763306214 0.271083990138 m 1 1 1 + +O #label +0 #magnetism +1 #number of atoms +0.517665952098 0.703037854903 0.321935083854 m 1 1 1 diff --git a/tests/abacus.relax/OUT.abacus/STRU_READIN_ADJUST.cif b/tests/abacus.relax/OUT.abacus/STRU_READIN_ADJUST.cif new file mode 100644 index 000000000..488d7efbd --- /dev/null +++ b/tests/abacus.relax/OUT.abacus/STRU_READIN_ADJUST.cif @@ -0,0 +1,22 @@ +data_test + +_audit_creation_method generated by ABACUS + +_cell_length_a 14.817 +_cell_length_b 14.817 +_cell_length_c 14.817 +_cell_angle_alpha 90 +_cell_angle_beta 90 +_cell_angle_gamma 90 + +_symmetry_space_group_name_H-M +_symmetry_Int_Tables_number + +loop_ +_atom_site_label +_atom_site_fract_x +_atom_site_fract_y +_atom_site_fract_z +H 0.569758 0.6702 0.29983 +H 0.491826 0.736268 0.271857 +O 0.518329 0.703007 0.31994 diff --git a/tests/abacus.relax/OUT.abacus/running_cell-relax.log b/tests/abacus.relax/OUT.abacus/running_cell-relax.log new file mode 100644 index 000000000..3fe32f482 --- /dev/null +++ b/tests/abacus.relax/OUT.abacus/running_cell-relax.log @@ -0,0 +1,1187 @@ + + WELCOME TO ABACUS + + 'Atomic-orbital Based Ab-initio Computation at UStc' + + Website: http://abacus.ustc.edu.cn/ + + Version: Parallel, in development + Processor Number is 2 + Start Time is Mon Jul 25 11:30:20 2022 + + ------------------------------------------------------------------------------------ + + READING GENERAL INFORMATION + global_out_dir = OUT.abacus/ + global_in_card = INPUT + pseudo_dir = + orbital_dir = + pseudo_type = auto + DRANK = 1 + DSIZE = 2 + DCOLOR = 1 + GRANK = 1 + GSIZE = 1 + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Reading atom information in unitcell: | + | From the input file and the structure file we know the number of | + | different elments in this unitcell, then we list the detail | + | information for each element, especially the zeta and polar atomic | + | orbital number for each element. The total atom number is counted. | + | We calculate the nearest atom distance for each atom and show the | + | Cartesian and Direct coordinates for each atom. We list the file | + | address for atomic orbitals. The volume and the lattice vectors | + | in real and reciprocal space is also shown. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + READING UNITCELL INFORMATION + ntype = 2 + atom label for species 1 = H + atom label for species 2 = O + lattice constant (Bohr) = 1 + lattice constant (Angstrom) = 0.529177 + + READING ATOM TYPE 1 + atom label = H + L=0, number of zeta = 1 + L=1, number of zeta = 1 + L=2, number of zeta = 1 + number of atom for this type = 2 + start magnetization = FALSE + start magnetization = FALSE + + READING ATOM TYPE 2 + atom label = O + L=0, number of zeta = 1 + L=1, number of zeta = 1 + L=2, number of zeta = 1 + number of atom for this type = 1 + start magnetization = FALSE + + TOTAL ATOM NUMBER = 3 + + CARTESIAN COORDINATES ( UNIT = 1 Bohr ). + atom x y z mag vx vy vz + tauc_H1 15.9532129411 18.7655861467 8.39524747132 0 0 0 0 + tauc_H2 13.7711312041 20.6154930027 7.61198952454 0 0 0 0 + tauc_O1 14.5132108826 19.6841922084 8.95832135273 0 0 0 0 + + + Volume (Bohr^3) = 21952 + Volume (A^3) = 3252.94689686 + + Lattice vectors: (Cartesian coordinate: in unit of a_0) + +28 +0 +0 + +0 +28 +0 + +0 +0 +28 + Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0) + +0.0357142857143 -0 +0 + +0 +0.0357142857143 -0 + +0 -0 +0.0357142857143 + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Reading pseudopotentials files: | + | The pseudopotential file is in UPF format. The 'NC' indicates that | + | the type of pseudopotential is 'norm conserving'. Functional of | + | exchange and correlation is decided by 4 given parameters in UPF | + | file. We also read in the 'core correction' if there exists. | + | Also we can read the valence electrons number and the maximal | + | angular momentum used in this pseudopotential. We also read in the | + | trail wave function, trail atomic density and local-pseudopotential| + | on logrithmic grid. The non-local pseudopotential projector is also| + | read in if there is any. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + PAO radial cut off (Bohr) = 15 + + Read in pseudopotential file is ../potential/H_ONCV_PBE-1.0.upf + pseudopotential type = NC + exchange-correlation functional = PBE + nonlocal core correction = 0 + valence electrons = 1 + lmax = 0 + number of zeta = 0 + number of projectors = 2 + L of projector = 0 + L of projector = 0 + PAO radial cut off (Bohr) = 15 + + Read in pseudopotential file is ../potential/O_ONCV_PBE-1.0.upf + pseudopotential type = NC + exchange-correlation functional = PBE + nonlocal core correction = 0 + valence electrons = 6 + lmax = 1 + number of zeta = 0 + number of projectors = 4 + L of projector = 0 + L of projector = 0 + L of projector = 1 + L of projector = 1 + initial pseudo atomic orbital number = 0 + NLOCAL = 27 + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Setup plane waves of charge/potential: | + | Use the energy cutoff and the lattice vectors to generate the | + | dimensions of FFT grid. The number of FFT grid on each processor | + | is 'nrxx'. The number of plane wave basis in reciprocal space is | + | different for charege/potential and wave functions. We also set | + | the 'sticks' for the parallel of FFT. The number of plane waves | + | is 'npw' in each processor. | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + + SETUP THE PLANE WAVE BASIS + energy cutoff for charge/potential (unit:Ry) = 200 + [fft grid for charge/potential] = 128, 128, 128 + [fft grid division] = 1, 1, 1 + [big fft grid for charge/potential] = 128, 128, 128 + nbxx = 1048576 + nrxx = 1048576 + + SETUP PLANE WAVES FOR CHARGE/POTENTIAL + number of plane waves = 1048171 + number of sticks = 12469 + + PARALLEL PW FOR CHARGE/POTENTIAL + PROC COLUMNS(POT) PW + 1 6235 524087 + 2 6234 524084 + --------------- sum ------------------- + 2 12469 1048171 + number of |g| = 3312 + max |g| = 5.06505102041 + min |g| = 0 + + SETUP THE ELECTRONS NUMBER + electron number of element H = 1 + total electron number of element H = 2 + electron number of element O = 6 + total electron number of element O = 6 + occupied bands = 4 + NBANDS = 6 + DONE : SETUP UNITCELL Time : 0.128108929377 (SEC) + + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Doing symmetry analysis: | + | We calculate the norm of 3 vectors and the angles between them, | + | the type of Bravais lattice is given. We can judge if the unticell | + | is a primitive cell. Finally we give the point group operation for | + | this unitcell. We we use the point group operations to do symmetry | + | analysis on given k-point mesh and the charge density. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + LATTICE VECTORS: (CARTESIAN COORDINATE: IN UNIT OF A0) + +28 +0 +0 + +0 +28 +0 + +0 +0 +28 + right hand lattice = 1 + NORM_A = 28 + NORM_B = 28 + NORM_C = 28 + ALPHA (DEGREE) = 90 + BETA (DEGREE) = 90 + GAMMA (DEGREE) = 90 + BRAVAIS TYPE = 1 + BRAVAIS LATTICE NAME = 01. Cubic P (simple) + IBRAV = 1 + BRAVAIS = SIMPLE CUBIC + LATTICE CONSTANT A = 122.049170419 + ibrav = 1 + ROTATION MATRICES = 48 + PURE POINT GROUP OPERATIONS = 1 + SPACE GROUP OPERATIONS = 1 + POINT GROUP = C_1 +Warning : If the optimal symmetric configuration is not the input configuration, +you have to manually change configurations, ABACUS would only calculate the input structure! + DONE : SYMMETRY Time : 0.155216244515 (SEC) + + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Setup K-points | + | We setup the k-points according to input parameters. | + | The reduced k-points are set according to symmetry operations. | + | We treat the spin as another set of k-points. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + + SETUP K-POINTS + nspin = 1 + Input type of k points = Monkhorst-Pack(Gamma) + nkstot = 1 + nkstot_ibz = 1 + IBZ DirectX DirectY DirectZ Weight ibz2bz + 1 0 0 0 1 0 + nkstot now = 1 + + KPOINTS DIRECT_X DIRECT_Y DIRECT_Z WEIGHT + 1 0 0 0 1 + + k-point number in this process = 1 + minimum distributed K point number = 1 + + KPOINTS CARTESIAN_X CARTESIAN_Y CARTESIAN_Z WEIGHT + 1 0 0 0 2 + + KPOINTS DIRECT_X DIRECT_Y DIRECT_Z WEIGHT + 1 0 0 0 2 + DONE : INIT K-POINTS Time : 0.155574926175 (SEC) + + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Setup plane waves of wave functions: | + | Use the energy cutoff and the lattice vectors to generate the | + | dimensions of FFT grid. The number of FFT grid on each processor | + | is 'nrxx'. The number of plane wave basis in reciprocal space is | + | different for charege/potential and wave functions. We also set | + | the 'sticks' for the parallel of FFT. The number of plane wave of | + | each k-point is 'npwk[ik]' in each processor | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + + SETUP PLANE WAVES FOR WAVE FUNCTIONS + energy cutoff for wavefunc (unit:Ry) = 50 + [fft grid for wave functions] = 128, 128, 128 + number of plane waves = 131155 + number of sticks = 3125 + + PARALLEL PW FOR WAVE FUNCTIONS + PROC COLUMNS(POT) PW + 1 1562 65576 + 2 1563 65579 + --------------- sum ------------------- + 2 3125 131155 + DONE : INIT PLANEWAVE Time : 0.174154018052 (SEC) + + DONE : INIT CHARGE Time : 0.259490167722 (SEC) + + npwx = 65576 + + SETUP NONLOCAL PSEUDOPOTENTIALS IN PLANE WAVE BASIS + H non-local projectors: + projector 1 L=0 + projector 2 L=0 + O non-local projectors: + projector 1 L=0 + projector 2 L=0 + projector 3 L=1 + projector 4 L=1 + TOTAL NUMBER OF NONLOCAL PROJECTORS = 12 + DONE : LOCAL POTENTIAL Time : 0.320586761925 (SEC) + + + Init Non-Local PseudoPotential table : + Init Non-Local-Pseudopotential done. + DONE : NON-LOCAL POTENTIAL Time : 0.331188020762 (SEC) + + init_chg = atomic + DONE : INIT POTENTIAL Time : 0.953154 (SEC) + + + Make real space PAO into reciprocal space. + max mesh points in Pseudopotential = 601 + dq(describe PAO in reciprocal space) = 0.01 + max q = 854 + + number of pseudo atomic orbitals for H is 0 + + number of pseudo atomic orbitals for O is 0 + DONE : INIT BASIS Time : 0.953302 (SEC) + + + ------------------------------------------- + RELAX CELL : 1 + RELAX IONS : 1 (in total: 1) + ------------------------------------------- + + PW ALGORITHM --------------- ION= 1 ELEC= 1-------------------------------- + + Density error is 0.417176464796 + Error Threshold = 0.01 + + Energy Rydberg eV + E_KohnSham -34.1543953066 -464.694387914 + E_Harris -34.3635778924 -467.540463003 + E_Fermi -0.520811172492 -7.08599952795 + + PW ALGORITHM --------------- ION= 1 ELEC= 2-------------------------------- + + Density error is 0.0243970602606 + Error Threshold = 0.00521470580995 + + Energy Rydberg eV + E_KohnSham -34.2350427506 -465.791652681 + E_Harris -34.2396178942 -465.853900704 + E_Fermi -0.445620982249 -6.06298450694 + + PW ALGORITHM --------------- ION= 1 ELEC= 3-------------------------------- + + Density error is 0.0106666914925 + Error Threshold = 0.000304963253257 + + Energy Rydberg eV + E_KohnSham -34.2334330467 -465.769751537 + E_Harris -34.2371193473 -465.819906229 + E_Fermi -0.461726538377 -6.28211183974 + + PW ALGORITHM --------------- ION= 1 ELEC= 4-------------------------------- + + Density error is 0.000502675170383 + Error Threshold = 0.000133333643656 + + Energy Rydberg eV + E_KohnSham -34.2339527078 -465.776821889 + E_Harris -34.2341613373 -465.779660439 + E_Fermi -0.200209037333 -2.72398369882 + + PW ALGORITHM --------------- ION= 1 ELEC= 5-------------------------------- + + Density error is 0.00013515778285 + Error Threshold = 6.28343962979e-06 + + Energy Rydberg eV + E_KohnSham -34.2339755848 -465.777133147 + E_Harris -34.2340600257 -465.778282024 + E_Fermi -0.1998556995 -2.71917629098 + + PW ALGORITHM --------------- ION= 1 ELEC= 6-------------------------------- + + Density error is 4.49530417282e-06 + Error Threshold = 1.68947228562e-06 + + Energy Rydberg eV + E_KohnSham -34.2340060994 -465.777548319 + E_Harris -34.2340025802 -465.777500437 + E_Fermi -0.198059680486 -2.69474019867 + + PW ALGORITHM --------------- ION= 1 ELEC= 7-------------------------------- + + Density error is 8.77501268413e-06 + Error Threshold = 5.61913021602e-08 + + Energy Rydberg eV + E_KohnSham -34.2340033243 -465.777510561 + E_Harris -34.2340088507 -465.777585752 + E_Fermi -0.197734942978 -2.69032191821 + + PW ALGORITHM --------------- ION= 1 ELEC= 8-------------------------------- + + Density error is 8.67497506757e-08 + Error Threshold = 5.61913021602e-08 + + Energy Rydberg eV + E_KohnSham -34.2340048296 -465.777531042 + E_Harris -34.234005046 -465.777533987 + E_band -8.10804363451 -110.315593062 + E_one_elec -69.5144554731 -945.792687802 + E_Hartree +36.17757886 +492.221212341 + E_xc -8.44392246259 -114.885458961 + E_Ewald +7.5467942461 +102.679403381 + E_demet -3.01469596081e-24 -4.10170428046e-23 + E_descf +0 +0 + E_efield +0 +0 + E_exx +0 +0 + E_Fermi -0.197510647395 -2.68727022024 + + charge density convergence is achieved + final etot is -465.777531042 eV + + STATE ENERGY(eV) AND OCCUPATIONS NSPIN == 1 + 1/1 kpoint (Cartesian) = 0 0 0 (65576 pws) + 1 -25.4315 2.00000 + 2 -13.5818 2.00000 + 3 -8.97266 2.00000 + 4 -7.17178 2.00000 + 5 -0.769414 0.00000 + 6 0.0954287 0.00000 + + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + H1 +0.40678204 -0.13991970 -0.61593726 + H2 +0.05888465 +0.16630779 -0.76223540 + O1 -0.46566669 -0.02638809 +1.37817266 + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-STRESS (KBAR) + + ><><><><><><><><><><><><><><><><><><><><><>< + + -2.153904 -0.330883 -0.070536 + -0.330883 -2.338728 -0.128010 + -0.070536 -0.128010 -1.978256 + TOTAL-PRESSURE: -2.156963 KBAR + + Ion relaxation is not converged yet (threshold is +0.025711) + + CARTESIAN COORDINATES ( UNIT = +1.000000 Bohr ). + atom x y z mag vx vy vz + tauc_H1 +16.065287158181 +18.727036287001 +8.225548029745 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + tauc_H2 +13.787354759370 +20.661313156171 +7.401982872047 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + tauc_O1 +14.384913110258 +19.676921914754 +9.338027446795 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + + Setup the structure factor in plane wave basis. + Setup the extrapolated charge. + NEW-OLD atomic charge density approx. for the potential ! + Setup the Vl+Vh+Vxc according to new structure factor and new charge. + Setup the new wave functions? + + ------------------------------------------- + RELAX CELL : 1 + RELAX IONS : 2 (in total: 2) + ------------------------------------------- + + PW ALGORITHM --------------- ION= +2 ELEC= +1-------------------------------- + + Density error is +4.306406951178 + Error Threshold = +0.010000000000 + + Energy Rydberg eV + E_KohnSham -33.984845752223 -462.387547881332 + E_Harris -36.146219524498 -491.794546692019 + E_Fermi -0.197995181756 -2.693862648432 + + PW ALGORITHM --------------- ION= 2 ELEC= 2-------------------------------- + + Density error is 0.524139822667 + Error Threshold = 0.010000000000 + + Energy Rydberg eV + E_KohnSham -34.028968257734 -462.987865366314 + E_Harris -34.358907885553 -467.476924300648 + E_Fermi -0.181612940365 -2.470970819494 + + PW ALGORITHM --------------- ION= 2 ELEC= 3-------------------------------- + + Density error is 0.008603629609 + Error Threshold = 0.006551747783 + + Energy Rydberg eV + E_KohnSham -34.126866845044 -464.319843979885 + E_Harris -34.126456524090 -464.314261276898 + E_Fermi -0.297352672177 -4.045690657129 + + PW ALGORITHM --------------- ION= 2 ELEC= 4-------------------------------- + + Density error is 0.002218320093 + Error Threshold = 0.000107545370 + + Energy Rydberg eV + E_KohnSham -34.129390206695 -464.354176076454 + E_Harris -34.129193594956 -464.351501036512 + E_Fermi -0.416104088025 -5.661386558240 + + PW ALGORITHM --------------- ION= 2 ELEC= 5-------------------------------- + + Density error is 0.001468898398 + Error Threshold = 0.000027729001 + + Energy Rydberg eV + E_KohnSham -34.129256085952 -464.352351270125 + E_Harris -34.130075919228 -464.363505674089 + E_Fermi -0.418170511223 -5.689501688208 + + PW ALGORITHM --------------- ION= 2 ELEC= 6-------------------------------- + + Density error is 0.000129882152 + Error Threshold = 0.000018361230 + + Energy Rydberg eV + E_KohnSham -34.129606602018 -464.357120285860 + E_Harris -34.129674949307 -464.358050198437 + E_Fermi -0.300798600818 -4.092574921552 + + PW ALGORITHM --------------- ION= 2 ELEC= 7-------------------------------- + + Density error is 0.000010859629 + Error Threshold = 0.000001623527 + + Energy Rydberg eV + E_KohnSham -34.129637781039 -464.357544498211 + E_Harris -34.129643192740 -464.357618128180 + E_Fermi -0.318223545585 -4.329653457719 + + PW ALGORITHM --------------- ION= 2 ELEC= 8-------------------------------- + + Density error is 0.000000556298 + Error Threshold = 0.000000135745 + + Energy Rydberg eV + E_KohnSham -34.129640373324 -464.357579768059 + E_Harris -34.129641983060 -464.357601669637 + E_Fermi -0.318727582976 -4.336511238241 + + PW ALGORITHM --------------- ION= 2 ELEC= 9-------------------------------- + + Density error is 0.000000121833 + Error Threshold = 0.000000006954 + + Energy Rydberg eV + E_KohnSham -34.129640005864 -464.357574768506 + E_Harris -34.129640436568 -464.357580628526 + E_Fermi -0.318641075078 -4.335334237908 + + PW ALGORITHM --------------- ION= 2 ELEC= 10-------------------------------- + + Density error is 0.000000186810 + Error Threshold = 0.000000001523 + + Energy Rydberg eV + E_KohnSham -34.129640196409 -464.357577360995 + E_Harris -34.129640113063 -464.357576227022 + E_Fermi -0.318550723465 -4.334104941148 + + PW ALGORITHM --------------- ION= 2 ELEC= 11-------------------------------- + + Density error is 0.000000001485 + Error Threshold = 0.000000001523 + + Energy Rydberg eV + E_KohnSham -34.129640061126 -464.357575520378 + E_Harris -34.129640228306 -464.357577794985 + E_band -7.520787696202 -102.325566116640 + E_one_elec -64.823282231884 -881.966001415785 + E_Hartree +33.884562763382 +461.023127820625 + E_xc -8.077485864562 -109.899833272499 + E_Ewald +4.886565271938 +66.485131347281 + E_demet -0.000000000000 -0.000000000000 + E_descf +0.000000000000 +0.000000000000 + E_efield +0.000000000000 +0.000000000000 + E_exx +0.000000000000 +0.000000000000 + E_Fermi -0.318592256747 -4.334670030437 + + charge density convergence is achieved + final etot is -464.357575520378 eV + + STATE ENERGY(eV) AND OCCUPATIONS NSPIN == 1 + 1/1 kpoint (Cartesian) = 0.00000 0.00000 0.00000 (65576 pws) + 1 -23.732691 2.000000 + 2 -10.717782 2.000000 + 3 -9.827251 2.000000 + 4 -6.885060 2.000000 + 5 -1.826406 0.000000 + 6 -0.276841 0.000000 + + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + H1 -3.43125048 +1.80114652 +2.82460187 + H2 +0.90733216 -1.89495640 +4.43921704 + O1 +2.52391832 +0.09380988 -7.26381891 + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-STRESS (KBAR) + + ><><><><><><><><><><><><><><><><><><><><><>< + + -4.009216 +0.908593 +0.503260 + +0.908593 -3.391991 +0.372926 + +0.503260 +0.372926 -5.317856 + TOTAL-PRESSURE: -4.239688 KBAR + + Ion relaxation is not converged yet (threshold is +0.025711) + + CARTESIAN COORDINATES ( UNIT = +1.000000 Bohr ). + atom x y z mag vx vy vz + tauc_H1 +15.970368922751 +18.759685051002 +8.369270394915 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + tauc_H2 +13.773614656790 +20.622507013014 +7.579842342394 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + tauc_O1 +14.493571448267 +19.683079293910 +9.016445611278 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + + Setup the structure factor in plane wave basis. + Setup the extrapolated charge. + NEW-OLD atomic charge density approx. for the potential ! + Setup the Vl+Vh+Vxc according to new structure factor and new charge. + Setup the new wave functions? + + ------------------------------------------- + RELAX CELL : 1 + RELAX IONS : 3 (in total: 3) + ------------------------------------------- + + PW ALGORITHM --------------- ION= +3 ELEC= +1-------------------------------- + + Density error is +2.863901623074 + Error Threshold = +0.010000000000 + + Energy Rydberg eV + E_KohnSham -34.163066102874 -464.812360149745 + E_Harris -35.594866834916 -484.293008506088 + E_Fermi -0.455725359897 -6.200461617704 + + PW ALGORITHM --------------- ION= 3 ELEC= 2-------------------------------- + + Density error is 0.410367630415 + Error Threshold = 0.010000000000 + + Energy Rydberg eV + E_KohnSham -34.147651136624 -464.602628774257 + E_Harris -34.409877454265 -468.170400859739 + E_Fermi -0.502269645612 -6.833729112759 + + PW ALGORITHM --------------- ION= 3 ELEC= 3-------------------------------- + + Density error is 0.016431897145 + Error Threshold = 0.005129595380 + + Energy Rydberg eV + E_KohnSham -34.231884315653 -465.748679969711 + E_Harris -34.237205015345 -465.821071802865 + E_Fermi -0.163677941108 -2.226952635976 + + PW ALGORITHM --------------- ION= 3 ELEC= 4-------------------------------- + + Density error is 0.003339586249 + Error Threshold = 0.000205398714 + + Energy Rydberg eV + E_KohnSham -34.234286271013 -465.781360248946 + E_Harris -34.234645115958 -465.786242584902 + E_Fermi -0.419282058198 -5.704625060666 + + PW ALGORITHM --------------- ION= 3 ELEC= 5-------------------------------- + + Density error is 0.007193242746 + Error Threshold = 0.000041744828 + + Energy Rydberg eV + E_KohnSham -34.233965842505 -465.777000595434 + E_Harris -34.237258750651 -465.821802909220 + E_Fermi -0.445575876230 -6.062370808069 + + PW ALGORITHM --------------- ION= 3 ELEC= 6-------------------------------- + + Density error is 0.000071123739 + Error Threshold = 0.000041744828 + + Energy Rydberg eV + E_KohnSham -34.235105030292 -465.792500040429 + E_Harris -34.235135970645 -465.792921005526 + E_Fermi -0.195917059478 -2.665588344306 + + PW ALGORITHM --------------- ION= 3 ELEC= 7-------------------------------- + + Density error is 0.000003025816 + Error Threshold = 0.000000889047 + + Energy Rydberg eV + E_KohnSham -34.235144953777 -465.793043227315 + E_Harris -34.235146413266 -465.793063084684 + E_Fermi -0.195891623194 -2.665242265913 + + PW ALGORITHM --------------- ION= 3 ELEC= 8-------------------------------- + + Density error is 0.000014047622 + Error Threshold = 0.000000037823 + + Energy Rydberg eV + E_KohnSham -34.235144028617 -465.793030639864 + E_Harris -34.235149064326 -465.793099154201 + E_Fermi -0.195695384814 -2.662572305769 + + PW ALGORITHM --------------- ION= 3 ELEC= 9-------------------------------- + + Density error is 0.000000380301 + Error Threshold = 0.000000037823 + + Energy Rydberg eV + E_KohnSham -34.235146948037 -465.793070360619 + E_Harris -34.235147173255 -465.793073424855 + E_Fermi -0.195885679726 -2.665161400881 + + PW ALGORITHM --------------- ION= 3 ELEC= 10-------------------------------- + + Density error is 0.000000003456 + Error Threshold = 0.000000004754 + + Energy Rydberg eV + E_KohnSham -34.235147175798 -465.793073459458 + E_Harris -34.235147030585 -465.793071483737 + E_band -8.005203857309 -108.916386110984 + E_one_elec -68.764358752319 -935.587098347715 + E_Hartree +35.812701320327 +487.256798728570 + E_xc -8.382702179786 -114.052514282117 + E_Ewald +7.099212435981 +96.589740441803 + E_demet -0.000000000000 -0.000000000000 + E_descf +0.000000000000 +0.000000000000 + E_efield +0.000000000000 +0.000000000000 + E_exx +0.000000000000 +0.000000000000 + E_Fermi -0.195843647331 -2.664589520805 + + charge density convergence is achieved + final etot is -465.793073459458 eV + + STATE ENERGY(eV) AND OCCUPATIONS NSPIN == 1 + 1/1 kpoint (Cartesian) = 0.00000 0.00000 0.00000 (65576 pws) + 1 -25.151687 2.000000 + 2 -13.061167 2.000000 + 3 -9.128130 2.000000 + 4 -7.117209 2.000000 + 5 -0.875259 0.000000 + 6 0.081372 0.000000 + + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + H1 -0.74519674 +0.51685204 +0.05757714 + H2 +0.52845057 -0.57175572 +0.56392399 + O1 +0.21674617 +0.05490368 -0.62150113 + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-STRESS (KBAR) + + ><><><><><><><><><><><><><><><><><><><><><>< + + -2.676771 +0.070442 -0.080953 + +0.070442 -2.674219 +0.037379 + -0.080953 +0.037379 -2.531715 + TOTAL-PRESSURE: -2.627568 KBAR + + Ion relaxation is not converged yet (threshold is +0.025711) + + CARTESIAN COORDINATES ( UNIT = +1.000000 Bohr ). + atom x y z mag vx vy vz + tauc_H1 +15.963906084978 +18.761908055509 +8.379056231391 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + tauc_H2 +13.772679114178 +20.619864761958 +7.591952522256 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + tauc_O1 +14.500969828653 +19.683498540459 +8.994549594941 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + + Setup the structure factor in plane wave basis. + Setup the extrapolated charge. + NEW-OLD atomic charge density approx. for the potential ! + Setup the Vl+Vh+Vxc according to new structure factor and new charge. + Setup the new wave functions? + + ------------------------------------------- + RELAX CELL : 1 + RELAX IONS : 4 (in total: 4) + ------------------------------------------- + + PW ALGORITHM --------------- ION= +4 ELEC= +1-------------------------------- + Notice: Threshold on eigenvalues was too large. + hsover_error=+0.080000 > DRHO=+0.009754 + Origin diag_ethr = +0.010000 + New diag_ethr = +0.000122 + + Density error is +0.013008603075 + Error Threshold = +0.000121928257 + + Energy Rydberg eV + E_KohnSham -34.235665065186 -465.800119706069 + E_Harris -34.242123145811 -465.887986400714 + E_Fermi -0.201851572452 -2.746331535609 + + PW ALGORITHM --------------- ION= 4 ELEC= 2-------------------------------- + + Density error is 0.002090692286 + Error Threshold = 0.000162607538 + + Energy Rydberg eV + E_KohnSham -34.235132150914 -465.792869035430 + E_Harris -34.236530970810 -465.811900956494 + E_Fermi -0.198476768045 -2.700414966033 + + PW ALGORITHM --------------- ION= 4 ELEC= 3-------------------------------- + + Density error is 0.000052460392 + Error Threshold = 0.000026133654 + + Energy Rydberg eV + E_KohnSham -34.235682052400 -465.800350828976 + E_Harris -34.235704411445 -465.800655039388 + E_Fermi -0.196153271944 -2.668802179785 + + PW ALGORITHM --------------- ION= 4 ELEC= 4-------------------------------- + + Density error is 0.000009658088 + Error Threshold = 0.000000655755 + + Energy Rydberg eV + E_KohnSham -34.235696222774 -465.800543626804 + E_Harris -34.235698393982 -465.800573167609 + E_Fermi -0.196616103423 -2.675099325104 + + PW ALGORITHM --------------- ION= 4 ELEC= 5-------------------------------- + + Density error is 0.000000187065 + Error Threshold = 0.000000120726 + + Energy Rydberg eV + E_KohnSham -34.235697732989 -465.800564174335 + E_Harris -34.235697255523 -465.800557678071 + E_Fermi -0.196471698719 -2.673134598311 + + PW ALGORITHM --------------- ION= 4 ELEC= 6-------------------------------- + + Density error is 0.000000075160 + Error Threshold = 0.000000002338 + + Energy Rydberg eV + E_KohnSham -34.235698022205 -465.800568109314 + E_Harris -34.235697872683 -465.800566074967 + E_band -8.042398038358 -109.422438905691 + E_one_elec -69.045489549154 -939.412079067952 + E_Hartree +35.948698866316 +489.107140268040 + E_xc -8.405372202565 -114.360955765699 + E_Ewald +7.266464863199 +98.865326456297 + E_demet -0.000000000000 -0.000000000000 + E_descf +0.000000000000 +0.000000000000 + E_efield +0.000000000000 +0.000000000000 + E_exx +0.000000000000 +0.000000000000 + E_Fermi -0.196459651857 -2.672970692356 + + charge density convergence is achieved + final etot is -465.800568109314 eV + + STATE ENERGY(eV) AND OCCUPATIONS NSPIN == 1 + 1/1 kpoint (Cartesian) = 0.00000 0.00000 0.00000 (65576 pws) + 1 -25.255180 2.000000 + 2 -13.252315 2.000000 + 3 -9.068312 2.000000 + 4 -7.135412 2.000000 + 5 -0.833001 0.000000 + 6 0.086578 0.000000 + + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + H1 -0.34425976 +0.29101536 -0.19310908 + H2 +0.37384222 -0.31991297 +0.09092313 + O1 -0.02958246 +0.02889761 +0.10218595 + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-STRESS (KBAR) + + ><><><><><><><><><><><><><><><><><><><><><>< + + -2.495529 -0.069759 -0.084163 + -0.069759 -2.561022 -0.020268 + -0.084163 -0.020268 -2.324147 + TOTAL-PRESSURE: -2.460232 KBAR + + Ion relaxation is not converged yet (threshold is +0.025711) + + CARTESIAN COORDINATES ( UNIT = +1.000000 Bohr ). + atom x y z mag vx vy vz + tauc_H1 +15.946191313792 +18.778838846668 +8.361024276793 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + tauc_H2 +13.796717055281 +20.601372573983 +7.590351723874 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + tauc_O1 +14.494646658736 +19.685059937275 +9.014182347920 +0.000000000000 +0.000000000000 +0.000000000000 +0.000000000000 + + Setup the structure factor in plane wave basis. + Setup the extrapolated charge. + NEW-OLD atomic charge density approx. for the potential ! + Setup the Vl+Vh+Vxc according to new structure factor and new charge. + Setup the new wave functions? + + ------------------------------------------- + RELAX CELL : 1 + RELAX IONS : 5 (in total: 5) + ------------------------------------------- + + PW ALGORITHM --------------- ION= +5 ELEC= +1-------------------------------- + Notice: Threshold on eigenvalues was too large. + hsover_error=+0.080000 > DRHO=+0.005280 + Origin diag_ethr = +0.010000 + New diag_ethr = +0.000066 + + Density error is +0.008353549675 + Error Threshold = +0.000065997708 + + Energy Rydberg eV + E_KohnSham -34.236235736125 -465.807884082530 + E_Harris -34.240337321619 -465.863689016079 + E_Fermi -0.438351093885 -5.964072601368 + + PW ALGORITHM --------------- ION= 5 ELEC= 2-------------------------------- + + Density error is 0.001090032164 + Error Threshold = 0.000104419371 + + Energy Rydberg eV + E_KohnSham -34.236307152675 -465.808855754541 + E_Harris -34.236930505462 -465.817336904304 + E_Fermi -0.196029016571 -2.667111598704 + + PW ALGORITHM --------------- ION= 5 ELEC= 3-------------------------------- + + Density error is 0.000082249260 + Error Threshold = 0.000013625402 + + Energy Rydberg eV + E_KohnSham -34.236527086358 -465.811848105802 + E_Harris -34.236560493475 -465.812302632946 + E_Fermi -0.196758447273 -2.677036012549 + + PW ALGORITHM --------------- ION= 5 ELEC= 4-------------------------------- + + Density error is 0.000130293705 + Error Threshold = 0.000001028116 + + Energy Rydberg eV + E_KohnSham -34.236562121478 -465.812324783072 + E_Harris -34.236605870085 -465.812920013399 + E_Fermi -0.197596470944 -2.688437909528 + + PW ALGORITHM --------------- ION= 5 ELEC= 5-------------------------------- + + Density error is 0.000024525955 + Error Threshold = 0.000001028116 + + Energy Rydberg eV + E_KohnSham -34.236557021747 -465.812255397668 + E_Harris -34.236571054853 -465.812446327870 + E_Fermi -0.197329109008 -2.684800263772 + + PW ALGORITHM --------------- ION= 5 ELEC= 6-------------------------------- + + Density error is 0.000001304759 + Error Threshold = 0.000000306574 + + Energy Rydberg eV + E_KohnSham -34.236563844352 -465.812348223977 + E_Harris -34.236563957567 -465.812349764348 + E_Fermi -0.197072028110 -2.681302498717 + + PW ALGORITHM --------------- ION= 5 ELEC= 7-------------------------------- + + Density error is 0.000000346081 + Error Threshold = 0.000000016309 + + Energy Rydberg eV + E_KohnSham -34.236564300744 -465.812354433499 + E_Harris -34.236564193956 -465.812352980576 + E_Fermi -0.197125544263 -2.682030623327 + + PW ALGORITHM --------------- ION= 5 ELEC= 8-------------------------------- + + Density error is 0.000000130954 + Error Threshold = 0.000000004326 + + Energy Rydberg eV + E_KohnSham -34.236564259962 -465.812353878642 + E_Harris -34.236564316674 -465.812354650242 + E_Fermi -0.197131965268 -2.682117985583 + + PW ALGORITHM --------------- ION= 5 ELEC= 9-------------------------------- + + Density error is 0.000000001123 + Error Threshold = 0.000000001637 + + Energy Rydberg eV + E_KohnSham -34.236564293143 -465.812354330090 + E_Harris -34.236564293813 -465.812354339209 + E_band -8.062071986763 -109.690116706159 + E_one_elec -69.117794401277 -940.395837049864 + E_Hartree +35.979373791465 +489.524494035791 + E_xc -8.410666291229 -114.432985537248 + E_Ewald +7.312522607898 +99.491974221231 + E_demet -0.000000000000 -0.000000000000 + E_descf +0.000000000000 +0.000000000000 + E_efield +0.000000000000 +0.000000000000 + E_exx +0.000000000000 +0.000000000000 + E_Fermi -0.197145167314 -2.682297608633 + + charge density convergence is achieved + final etot is -465.812354330090 eV + + STATE ENERGY(eV) AND OCCUPATIONS NSPIN == 1 + 1/1 kpoint (Cartesian) = 0.00000 0.00000 0.00000 (65576 pws) + 1 -25.330615 2.000000 + 2 -13.155666 2.000000 + 3 -9.203788 2.000000 + 4 -7.154990 2.000000 + 5 -0.831095 0.000000 + 6 0.086695 0.000000 + + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + H1 -0.18807867 +0.13569889 -0.06941882 + H2 +0.14926844 -0.14376688 +0.04694486 + O1 +0.03881023 +0.00806799 +0.02247396 + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-STRESS (KBAR) + + ><><><><><><><><><><><><><><><><><><><><><>< + + -2.392469 -0.160549 -0.032796 + -0.160549 -2.480798 -0.059783 + -0.032796 -0.059783 -2.315957 + TOTAL-PRESSURE: -2.396408 KBAR + + + -------------------------------------------- + !FINAL_ETOT_IS -465.8123543300895335 eV + -------------------------------------------- + + + + + + + |CLASS_NAME---------|NAME---------------|TIME(Sec)-----|CALLS----|AVG------|PER%------- + total +84.59445 19 +4.45 +100.00% + Run_pw plane_wave_line +84.58434 1 +84.58 +99.99% + PW_Basis setup_struc_factor +0.37204 9 +0.04 +0.44% + Potential init_pot +2.60817 5 +0.52 +3.08% + Potential set_local_pot +0.17850 5 +0.04 +0.21% + PW_Basis recip2real +6.88976 294 +0.02 +8.14% + PW_Basis gathers_scatterp +2.80189 294 +0.01 +3.31% + Charge atomic_rho +0.93431 9 +0.10 +1.10% + Potential v_of_rho +21.19138 49 +0.43 +25.05% + XC_Functional v_xc +19.67060 54 +0.36 +23.25% + PW_Basis real2recip +9.16010 446 +0.02 +10.83% + PW_Basis gatherp_scatters +3.68783 446 +0.01 +4.36% + H_Hartree_pw v_hartree +3.16563 49 +0.06 +3.74% + Potential set_vr_eff +0.12252 49 +0.00 +0.14% + Cell_PW opt_cells_pw +83.63640 1 +83.64 +98.87% + Ions opt_ions_pw +83.63639 1 +83.64 +98.87% + ESolver_KS_PW Run +73.96240 5 +14.79 +87.43% + Symmetry rho_symmetry +3.36747 51 +0.07 +3.98% + HSolverPW solve +43.08348 46 +0.94 +50.93% + pp_cell_vnl getvnl +1.24764 56 +0.02 +1.47% + WF_igk get_sk +0.31200 231 +0.00 +0.37% + DiagoIterAssist diagH_subspace +6.55659 45 +0.15 +7.75% + HamiltPW h_psi +34.80823 1083 +0.03 +41.15% + Operator EkineticPW +0.20066 1083 +0.00 +0.24% + Operator VeffPW +31.49160 1083 +0.03 +37.23% + PW_Basis_K recip2real +17.21913 1493 +0.01 +20.35% + PW_Basis_K gathers_scatterp +4.88147 1493 +0.00 +5.77% + PW_Basis_K real2recip +11.93248 1308 +0.01 +14.11% + PW_Basis_K gatherp_scatters +2.55499 1308 +0.00 +3.02% + Operator NonlocalPW +3.11107 1083 +0.00 +3.68% + NonlocalPW add_nonlocal_pp +1.47037 1083 +0.00 +1.74% + DiagoCG diag_once +31.41731 46 +0.68 +37.14% + ElecStatePW psiToRho +3.97135 46 +0.09 +4.69% + Charge rho_mpi +1.34662 46 +0.03 +1.59% + Charge mix_rho +3.23741 39 +0.08 +3.83% + Forces cal_force_loc +0.29882 5 +0.06 +0.35% + Forces cal_force_ew +0.22785 5 +0.05 +0.27% + Forces cal_force_nl +0.20171 5 +0.04 +0.24% + Stress_PW cal_stress +3.22335 5 +0.64 +3.81% + Stress_Func stress_har +0.17289 5 +0.03 +0.20% + Stress_Func stress_ew +0.28902 5 +0.06 +0.34% + Stress_Func stress_gga +0.80316 5 +0.16 +0.95% + Stress_Func stress_loc +0.48018 5 +0.10 +0.57% + Stress_Func stres_nl +1.42753 5 +0.29 +1.69% + ---------------------------------------------------------------------------------------- + + CLASS_NAME---------|NAME---------------|MEMORY(MB)-------- + +418.2583 + Charge_Pulay Rrho +64.0000 + Charge_Pulay dRrho +56.0000 + Charge_Pulay drho +56.0000 + PW_Basis struc_fac +15.9939 + Charge rho +8.0000 + Charge rho_save +8.0000 + Charge rho_core +8.0000 + Potential vltot +8.0000 + Potential vr +8.0000 + Potential vr_eff +8.0000 + Potential vr_eff1 +8.0000 + Potential vnew +8.0000 + Charge_Pulay rho_save2 +8.0000 + wavefunc psi +6.0037 + Charge rhog +3.9985 + Charge rhog_save +3.9985 + Charge kin_r +3.9985 + Charge kin_r_save +3.9985 + Charge rhog_core +3.9985 + ---------------------------------------------------------- + + Start Time : Mon Jul 25 11:30:20 2022 + Finish Time : Mon Jul 25 11:31:45 2022 + Total Time : 0 h 1 mins 25 secs diff --git a/tests/abacus.relax/STRU b/tests/abacus.relax/STRU new file mode 100644 index 000000000..c1ae18bf0 --- /dev/null +++ b/tests/abacus.relax/STRU @@ -0,0 +1,29 @@ +ATOMIC_SPECIES +H 1.008 ../potential/H_ONCV_PBE-1.0.upf +O 15.9994 ../potential/O_ONCV_PBE-1.0.upf + +NUMERICAL_ORBITAL +../potential/H_gga_8au_100Ry_2s1p.orb +../potential/O_gga_7au_100Ry_2s2p1d.orb + + +LATTICE_CONSTANT +1 + +LATTICE_VECTORS +28 0 0 +0 28 0 +0 0 28 + +ATOMIC_POSITIONS +Cartesian + +H +0 +2 +-12.046787058887078 18.76558614676448 8.395247471328744 1 1 1 +-14.228868795885418 20.61549300274637 7.611989524516571 1 1 1 +O +0 +1 +-13.486789117423204 19.684192208418636 8.958321352749174 1 1 1 diff --git a/tests/abacus.relax/coord.ref b/tests/abacus.relax/coord.ref new file mode 100644 index 000000000..9181ed3f6 --- /dev/null +++ b/tests/abacus.relax/coord.ref @@ -0,0 +1,15 @@ + 8.44207672911 9.93032053807 4.44257364171 + 7.28736880156 10.90924908856 4.02809138602 + 7.68006045610 10.41642593172 4.74053950781 + 8.50138385072 9.90992083083 4.35277256453 + 7.29595393729 10.93349606958 3.91696065138 + 7.61216819877 10.41257865801 4.94147131963 + 8.45115528363 9.92719781271 4.42882716487 + 7.28868298813 10.91296074297 4.01107982983 + 7.66966771502 10.41583700273 4.77129754083 + 8.44773529717 9.92837417603 4.43400560653 + 7.28818792030 10.91156252393 4.01748826104 + 7.67358276932 10.41605885845 4.75971066798 + 8.43836104396 9.93733356488 4.42446350709 + 7.30090825093 10.90177687947 4.01664115501 + 7.67023669189 10.41688511407 4.77009987344 diff --git a/tests/abacus.relax/force.ref b/tests/abacus.relax/force.ref new file mode 100644 index 000000000..1b68d00be --- /dev/null +++ b/tests/abacus.relax/force.ref @@ -0,0 +1,15 @@ + 0.40678204000 -0.13991970000 -0.61593726000 + 0.05888465000 0.16630779000 -0.76223540000 + -0.46566669000 -0.02638809000 1.37817266000 + -3.43125048000 1.80114652000 2.82460187000 + 0.90733216000 -1.89495640000 4.43921704000 + 2.52391832000 0.09380988000 -7.26381891000 + -0.74519674000 0.51685204000 0.05757714000 + 0.52845057000 -0.57175572000 0.56392399000 + 0.21674617000 0.05490368000 -0.62150113000 + -0.34425976000 0.29101536000 -0.19310908000 + 0.37384222000 -0.31991297000 0.09092313000 + -0.02958246000 0.02889761000 0.10218595000 + -0.18807867000 0.13569889000 -0.06941882000 + 0.14926844000 -0.14376688000 0.04694486000 + 0.03881023000 0.00806799000 0.02247396000 diff --git a/tests/abacus.relax/stress.ref b/tests/abacus.relax/stress.ref new file mode 100644 index 000000000..7ba6e6954 --- /dev/null +++ b/tests/abacus.relax/stress.ref @@ -0,0 +1,15 @@ + -2.15390400000 -0.33088300000 -0.07053600000 + -0.33088300000 -2.33872800000 -0.12801000000 + -0.07053600000 -0.12801000000 -1.97825600000 + -4.00921600000 0.90859300000 0.50326000000 + 0.90859300000 -3.39199100000 0.37292600000 + 0.50326000000 0.37292600000 -5.31785600000 + -2.67677100000 0.07044200000 -0.08095300000 + 0.07044200000 -2.67421900000 0.03737900000 + -0.08095300000 0.03737900000 -2.53171500000 + -2.49552900000 -0.06975900000 -0.08416300000 + -0.06975900000 -2.56102200000 -0.02026800000 + -0.08416300000 -0.02026800000 -2.32414700000 + -2.39246900000 -0.16054900000 -0.03279600000 + -0.16054900000 -2.48079800000 -0.05978300000 + -0.03279600000 -0.05978300000 -2.31595700000 diff --git a/tests/abacus.relax/virial.ref b/tests/abacus.relax/virial.ref new file mode 100644 index 000000000..cb442e805 --- /dev/null +++ b/tests/abacus.relax/virial.ref @@ -0,0 +1,15 @@ + -4.37314061483 -0.67180240440 -0.14321151101 + -0.67180240440 -4.74839473061 -0.25990282302 + -0.14321151101 -0.25990282302 -4.01651682718 + -8.14004028185 1.84474561106 1.02178497548 + 1.84474561106 -6.88686849889 0.75716366046 + 1.02178497548 0.75716366046 -10.79701419257 + -5.43473431346 0.14302065978 -0.16436148138 + 0.14302065978 -5.42955290573 0.07589178675 + -0.16436148138 0.07589178675 -5.14022244802 + -5.06675284757 -0.14163394290 -0.17087884770 + -0.14163394290 -5.19972539337 -0.04115077273 + -0.17087884770 -0.04115077273 -4.71879045702 + -4.85750681257 -0.32596780199 -0.06658677434 + -0.32596780199 -5.03684402415 -0.12137934902 + -0.06658677434 -0.12137934902 -4.70216203642 diff --git a/tests/test_abacus_relax.py b/tests/test_abacus_relax.py new file mode 100644 index 000000000..3f1230ad6 --- /dev/null +++ b/tests/test_abacus_relax.py @@ -0,0 +1,75 @@ +import os +import numpy as np +import unittest +from context import dpdata +from dpdata.unit import LengthConversion + +bohr2ang = LengthConversion("bohr", "angstrom").value() + +class TestABACUSRelax: + + def test_atom_names(self) : + self.assertEqual(self.system.data['atom_names'], ['H','O']) + + def test_atom_numbs(self) : + self.assertEqual(self.system.data['atom_numbs'], [2,1]) + + def test_atom_types(self) : + ref_type = np.array([0,0,1]) + for ii in range(ref_type.shape[0]) : + self.assertEqual(self.system.data['atom_types'][ii], ref_type[ii]) + + def test_cell(self) : + cell = bohr2ang * 28.0 * np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + for idx in range(np.shape(self.system.data['cells'])[0]): + np.testing.assert_almost_equal(cell, self.system.data['cells'][idx], decimal = 5) + + def test_coord(self) : + with open('abacus.relax/coord.ref') as fp: + ref = [] + for ii in fp : + ref.append([float(jj) for jj in ii.split()]) + ref = np.array(ref) + ref = ref.reshape([5, 3, 3]) + np.testing.assert_almost_equal(self.system.data['coords'], ref, decimal = 5) + + def test_force(self) : + with open('abacus.relax/force.ref') as fp: + ref = [] + for ii in fp : + ref.append([float(jj) for jj in ii.split()]) + ref = np.array(ref) + ref = ref.reshape([5, 3, 3]) + np.testing.assert_almost_equal(self.system.data['forces'], ref, decimal=5) + + def test_virial(self) : + with open('abacus.relax/virial.ref') as fp: + ref = [] + for ii in fp : + ref.append([float(jj) for jj in ii.split()]) + ref = np.array(ref) + ref = ref.reshape([5, 3, 3]) + np.testing.assert_almost_equal(self.system.data['virials'], ref, decimal=5) + + def test_stress(self) : + with open('abacus.relax/stress.ref') as fp: + ref = [] + for ii in fp : + ref.append([float(jj) for jj in ii.split()]) + ref = np.array(ref) + ref = ref.reshape([5, 3, 3]) + np.testing.assert_almost_equal(self.system.data['stress'], ref, decimal=5) + + def test_energy(self) : + ref_energy = np.array([-465.77753104, -464.35757552, -465.79307346, -465.80056811, + -465.81235433]) + np.testing.assert_almost_equal(self.system.data['energies'], ref_energy) + + +class TestABACUSMDLabeledOutput(unittest.TestCase, TestABACUSRelax): + + def setUp(self): + self.system = dpdata.LabeledSystem('abacus.relax',fmt='abacus/relax') + +if __name__ == '__main__': + unittest.main() \ No newline at end of file