From e12271d3647b9ec8e827990378bf5aadf0444a8e Mon Sep 17 00:00:00 2001 From: LiuRenxi Date: Wed, 8 Jun 2022 22:23:50 +0800 Subject: [PATCH 1/4] adapt abacus/md interface to MD output without stress calculation. --- dpdata/abacus/md.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/dpdata/abacus/md.py b/dpdata/abacus/md.py index 8bd9c72cf..32774ddca 100644 --- a/dpdata/abacus/md.py +++ b/dpdata/abacus/md.py @@ -33,8 +33,17 @@ def get_coord_dump_freq(inlines): def get_coords_from_dump(dumplines, natoms): nlines = len(dumplines) total_natoms = sum(natoms) - nframes_dump = int(nlines/(total_natoms + 13)) - + calc_stress = False + if "VIRIAL" in dumplines[6]: + calc_stress = True + else: + assert("POSITIONS" in dumplines[6] and "FORCE" in dumplines[6]) + nframes_dump = -1 + if calc_stress: + nframes_dump = int(nlines/(total_natoms + 13)) + else: + nframes_dump = int(nlines/(total_natoms + 9)) + assert(nframes_dump > 0) cells = np.zeros([nframes_dump, 3, 3]) stresses = np.zeros([nframes_dump, 3, 3]) forces = np.zeros([nframes_dump, total_natoms, 3]) @@ -47,10 +56,15 @@ def get_coords_from_dump(dumplines, natoms): # read in LATTICE_VECTORS for ix in range(3): cells[iframe, ix] = np.array([float(i) for i in re.split('\s+', dumplines[iline+3+ix])[-3:]]) * celldm - stresses[iframe, ix] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+ix])[-3:]]) + if calc_stress: + stresses[iframe, ix] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+ix])[-3:]]) for iat in range(total_natoms): - coords[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+11+iat])[-6:-3]])*celldm - forces[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+11+iat])[-3:]]) + if calc_stress: + coords[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+11+iat])[-6:-3]])*celldm + forces[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+11+iat])[-3:]]) + else: + coords[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+iat])[-6:-3]])*celldm + forces[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+iat])[-3:]]) iframe += 1 assert(iframe == nframes_dump) cells *= bohr2ang From 61e1040e09632766720789c0b5abf29f495313c3 Mon Sep 17 00:00:00 2001 From: LiuRenxi Date: Thu, 9 Jun 2022 13:46:16 +0800 Subject: [PATCH 2/4] add unit test for abacus md without stress output. --- tests/abacus.md.nostress/INPUT | 27 ++++++++++++ tests/abacus.md.nostress/OUT.autotest/MD_dump | 44 +++++++++++++++++++ tests/abacus.md.nostress/STRU | 19 ++++++++ tests/abacus.md.nostress/Si_coord | 8 ++++ tests/abacus.md.nostress/Si_force | 8 ++++ tests/test_abacus_md.py | 41 ++++++++++++++++- 6 files changed, 145 insertions(+), 2 deletions(-) create mode 100644 tests/abacus.md.nostress/INPUT create mode 100644 tests/abacus.md.nostress/OUT.autotest/MD_dump create mode 100644 tests/abacus.md.nostress/STRU create mode 100644 tests/abacus.md.nostress/Si_coord create mode 100644 tests/abacus.md.nostress/Si_force diff --git a/tests/abacus.md.nostress/INPUT b/tests/abacus.md.nostress/INPUT new file mode 100644 index 000000000..4dd5c0fa6 --- /dev/null +++ b/tests/abacus.md.nostress/INPUT @@ -0,0 +1,27 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ./ +ntype 1 +nbands 8 +calculation md + +#Parameters (Accuracy) +ecutwfc 20 +scf_nmax 20 + +basis_type pw +md_nstep 3 + +#cal_stress 1 +stress_thr 1e-6 +cal_force 1 +force_thr_ev 1.0e-3 + +ks_solver cg +mixing_type pulay +mixing_beta 0.7 + +md_type -1 +md_dt 1 +md_tfirst 0 diff --git a/tests/abacus.md.nostress/OUT.autotest/MD_dump b/tests/abacus.md.nostress/OUT.autotest/MD_dump new file mode 100644 index 000000000..12bde6abd --- /dev/null +++ b/tests/abacus.md.nostress/OUT.autotest/MD_dump @@ -0,0 +1,44 @@ +MDSTEP: 0 +LATTICE_CONSTANT: 10.200000000000 +LATTICE_VECTORS + 0.500000000000 0.500000000000 0.000000000000 + 0.500000000000 0.000000000000 0.500000000000 + 0.000000000000 0.500000000000 0.500000000000 +INDEX LABEL POSITIONS FORCE (eV/Angstrom) + 0 Si 0.000000000000 0.000000000000 0.000000000000 -0.885362725259 0.500467424429 0.150239620221 + 1 Si 0.241000000000 0.255000000000 0.250999999999 0.885362725259 -0.500467424429 -0.150239620221 + + +MDSTEP: 1 +LATTICE_CONSTANT: 10.200000000000 +LATTICE_VECTORS + 0.500000000000 0.500000000000 0.000000000000 + 0.500000000000 0.000000000000 0.500000000000 + 0.000000000000 0.500000000000 0.500000000000 +INDEX LABEL POSITIONS FORCE (eV/Angstrom) + 0 Si 0.999208682339 0.500447306737 0.500134280856 -0.728804779648 0.408578746723 0.107042476463 + 1 Si 0.241791317661 0.254552693263 0.250865719143 0.728804779648 -0.408578746723 -0.107042476463 + + +MDSTEP: 2 +LATTICE_CONSTANT: 10.200000000000 +LATTICE_VECTORS + 0.500000000000 0.500000000000 0.000000000000 + 0.500000000000 0.000000000000 0.500000000000 + 0.000000000000 0.500000000000 0.500000000000 +INDEX LABEL POSITIONS FORCE (eV/Angstrom) + 0 Si 0.997114226602 0.501624803733 0.500458153134 -0.316038867402 0.171613491764 0.014802659803 + 1 Si 0.243885773398 0.253375196267 0.250541846865 0.316038867402 -0.171613491764 -0.014802659803 + + +MDSTEP: 3 +LATTICE_CONSTANT: 10.200000000000 +LATTICE_VECTORS + 0.500000000000 0.500000000000 0.000000000000 + 0.500000000000 0.000000000000 0.500000000000 + 0.000000000000 0.500000000000 0.500000000000 +INDEX LABEL POSITIONS FORCE (eV/Angstrom) + 0 Si 0.994451593844 0.503106810892 0.500786060548 0.204324467775 -0.117398116295 -0.052955519257 + 1 Si 0.246548406156 0.251893189108 0.250213939451 -0.204324467775 0.117398116295 0.052955519257 + + diff --git a/tests/abacus.md.nostress/STRU b/tests/abacus.md.nostress/STRU new file mode 100644 index 000000000..1533a93a6 --- /dev/null +++ b/tests/abacus.md.nostress/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Si 1 ../tools/PP_ORB/Si_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +10.2 + +LATTICE_VECTORS +0.5 0.5 0 #latvec1 +0.5 0 0.5 #latvec2 +0 0.5 0.5 #latvec3 + +ATOMIC_POSITIONS +Cartesian + +Si #label +0 #magnetism +2 #number of atoms +0 0 0 m 1 1 1 v -0 0 0 +0.241 0.255 0.251 m 1 1 1 v 0 -0 -0 diff --git a/tests/abacus.md.nostress/Si_coord b/tests/abacus.md.nostress/Si_coord new file mode 100644 index 000000000..177fd3690 --- /dev/null +++ b/tests/abacus.md.nostress/Si_coord @@ -0,0 +1,8 @@ +0 0 0 +1.30082342 1.37638993 1.3547995 +5.39333633 2.70121816 2.69952857 +1.30509464 1.37397554 1.3540747 +5.38203128 2.70757383 2.70127671 +1.31639969 1.36761987 1.35232656 +5.36765943 2.71557312 2.70304662 +1.33077154 1.35962058 1.35055665 \ No newline at end of file diff --git a/tests/abacus.md.nostress/Si_force b/tests/abacus.md.nostress/Si_force new file mode 100644 index 000000000..ab4566813 --- /dev/null +++ b/tests/abacus.md.nostress/Si_force @@ -0,0 +1,8 @@ +-0.88536273 0.50046742 0.15023962 +0.88536273 -0.50046742 -0.15023962 +-0.72880478 0.40857875 0.10704248 +0.72880478 -0.40857875 -0.10704248 +-0.31603887 0.17161349 0.01480266 +0.31603887 -0.17161349 -0.01480266 +0.20432447 -0.11739812 -0.05295552 +-0.20432447 0.11739812 0.05295552 \ No newline at end of file diff --git a/tests/test_abacus_md.py b/tests/test_abacus_md.py index 2ed61a994..2b40b3840 100644 --- a/tests/test_abacus_md.py +++ b/tests/test_abacus_md.py @@ -10,22 +10,32 @@ class TestABACUSMD: def test_atom_names(self) : self.assertEqual(self.system_water.data['atom_names'], ['H', 'O']) + self.assertEqual(self.system_Si.data['atom_names'], ['Si']) def test_atom_numbs(self) : self.assertEqual(self.system_water.data['atom_numbs'], [2, 1]) + self.assertEqual(self.system_Si.data['atom_numbs'], [2]) def test_atom_types(self) : ref_type = [0, 0, 1] ref_type = np.array(ref_type) + ref_type2 = np.array([0, 0]) for ii in range(ref_type.shape[0]) : self.assertEqual(self.system_water.data['atom_types'][ii], ref_type[ii]) + for ii in range(ref_type2.shape[0]) : + self.assertEqual(self.system_Si.data['atom_types'][ii], ref_type2[ii]) def test_cell(self) : cell = bohr2ang * 28 * np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + cell2 = bohr2ang * 5.1 * np.array([[1, 1, 0], [1, 0, 1], [0, 1, 1]]) for idx in range(np.shape(self.system_water.data['cells'])[0]): for ii in range(cell.shape[0]) : for jj in range(cell.shape[1]) : self.assertAlmostEqual(self.system_water.data['cells'][idx][ii][jj], cell[ii][jj]) + for idx in range(np.shape(self.system_Si.data['cells'])[0]): + for ii in range(cell2.shape[0]) : + for jj in range(cell2.shape[1]) : + self.assertAlmostEqual(self.system_Si.data['cells'][idx][ii][jj], cell2[ii][jj]) def test_coord(self) : fp = open('abacus.md/water_coord') @@ -40,6 +50,18 @@ def test_coord(self) : self.assertAlmostEqual(self.system_water.data['coords'][ii][jj][kk], coord[ii][jj][kk]) fp.close() + fp2 = open('abacus.md.nostress/Si_coord') + coord = [] + for ii in fp2 : + coord.append([float(jj) for jj in ii.split()]) + coord = np.array(coord) + coord = coord.reshape([4, 2, 3]) + for ii in range(coord.shape[0]) : + for jj in range(coord.shape[1]) : + for kk in range(coord.shape[2]): + self.assertAlmostEqual(self.system_Si.data['coords'][ii][jj][kk], coord[ii][jj][kk]) + fp2.close() + def test_force(self) : fp = open('abacus.md/water_force') force = [] @@ -53,6 +75,18 @@ def test_force(self) : self.assertAlmostEqual(self.system_water.data['forces'][ii][jj][kk], force[ii][jj][kk]) fp.close() + fp2 = open('abacus.md.nostress/Si_force') + force = [] + for ii in fp2 : + force.append([float(jj) for jj in ii.split()]) + force = np.array(force) + force = force.reshape([4, 2, 3]) + for ii in range(force.shape[0]) : + for jj in range(force.shape[1]) : + for kk in range(force.shape[2]): + self.assertAlmostEqual(self.system_Si.data['forces'][ii][jj][kk], force[ii][jj][kk]) + fp2.close() + def test_virial(self) : fp = open('abacus.md/water_virial') virial = [] @@ -69,15 +103,18 @@ def test_virial(self) : def test_energy(self) : ref_energy = np.array([-466.69285117, -466.69929051, -466.69829826, -466.70364664, -466.6976083]) + ref_energy2 = np.array([-211.77184603, -211.78111966, -211.79681663, -211.79875524]) for ii in range(5): self.assertAlmostEqual(self.system_water.data['energies'][ii], ref_energy[ii]) + for ii in range(4): + self.assertAlmostEqual(self.system_Si.data['energies'][ii], ref_energy2[ii]) class TestABACUSMDLabeledOutput(unittest.TestCase, TestABACUSMD): def setUp(self): - self.system_water = dpdata.LabeledSystem('abacus.md',fmt='abacus/md') - + self.system_water = dpdata.LabeledSystem('abacus.md',fmt='abacus/md') # system with stress + self.system_Si = dpdata.LabeledSystem('abacus.md.nostress',fmt='abacus/md') # system without stress if __name__ == '__main__': unittest.main() \ No newline at end of file From a5e42a7ddbafee560ccda80bde26008cfdfb20a3 Mon Sep 17 00:00:00 2001 From: Liu Renxi <75369672+Liu-RX@users.noreply.github.com> Date: Thu, 9 Jun 2022 13:51:05 +0800 Subject: [PATCH 3/4] Add files via upload --- .../OUT.autotest/running_md.log | 527 ++++++++++++++++++ 1 file changed, 527 insertions(+) create mode 100644 tests/abacus.md.nostress/OUT.autotest/running_md.log diff --git a/tests/abacus.md.nostress/OUT.autotest/running_md.log b/tests/abacus.md.nostress/OUT.autotest/running_md.log new file mode 100644 index 000000000..1f5269321 --- /dev/null +++ b/tests/abacus.md.nostress/OUT.autotest/running_md.log @@ -0,0 +1,527 @@ + + WELCOME TO ABACUS + + 'Atomic-orbital Based Ab-initio Computation at UStc' + + Website: http://abacus.ustc.edu.cn/ + + Version: Parallel, in development + Processor Number is 1 + Start Time is Thu Jun 9 13:13:06 2022 + + ------------------------------------------------------------------------------------ + + READING GENERAL INFORMATION + global_out_dir = OUT.autotest/ + global_in_card = INPUT + pseudo_dir = + orbital_dir = + pseudo_type = auto + DRANK = 1 + DSIZE = 1 + 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 = 1 + atom label for species 1 = Si + lattice constant (Bohr) = 10.2 + lattice constant (Angstrom) = 5.39761 + + READING ATOM TYPE 1 + atom label = Si + 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 + + TOTAL ATOM NUMBER = 2 + + CARTESIAN COORDINATES ( UNIT = 10.2 Bohr ). + atom x y z mag vx vy vz + tauc_Si1 0 0 0 0 -0 0 0 + tauc_Si2 0.241 0.255 0.250999999999 0 0 -0 -0 + + + Volume (Bohr^3) = 265.302 + Volume (A^3) = 39.3136533177 + + Lattice vectors: (Cartesian coordinate: in unit of a_0) + +0.5 +0.5 +0 + +0.5 +0 +0.5 + +0 +0.5 +0.5 + Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0) + +1 +1 -1 + +1 -1 +1 + -1 +1 +1 + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | 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 ../tools/PP_ORB/Si_ONCV_PBE-1.0.upf + pseudopotential type = NC + exchange-correlation functional = PBE + nonlocal core correction = 0 + valence electrons = 4 + 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 = 18 + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | 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) = 80 + use input fft dimensions for wave functions. + [fft grid for charge/potential] = 24, 24, 24 + [fft grid division] = 1, 1, 1 + [big fft grid for charge/potential] = 24, 24, 24 + nbxx = 13824 + nrxx = 13824 + + SETUP PLANE WAVES FOR CHARGE/POTENTIAL + number of plane waves = 3143 + number of sticks = 283 + + PARALLEL PW FOR CHARGE/POTENTIAL + PROC COLUMNS(POT) PW + 1 283 3143 + --------------- sum ------------------- + 1 283 3143 + number of |g| = 72 + max |g| = 208 + min |g| = 0 + + SETUP THE ELECTRONS NUMBER + electron number of element Si = 4 + total electron number of element Si = 8 + occupied bands = 4 + NBANDS = 8 + DONE : SETUP UNITCELL Time : 0.0560910534114 (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 = 8 + + KPOINTS DIRECT_X DIRECT_Y DIRECT_Z WEIGHT + 1 0 0 0 0.125 + 2 0.5 0 0 0.125 + 3 0 0.5 0 0.125 + 4 0.5 0.5 0 0.125 + 5 0 0 0.5 0.125 + 6 0.5 0 0.5 0.125 + 7 0 0.5 0.5 0.125 + 8 0.5 0.5 0.5 0.125 + + k-point number in this process = 8 + minimum distributed K point number = 8 + + KPOINTS CARTESIAN_X CARTESIAN_Y CARTESIAN_Z WEIGHT + 1 0 0 0 0.25 + 2 0.5 0.5 -0.5 0.25 + 3 0.5 -0.5 0.5 0.25 + 4 1 0 0 0.25 + 5 -0.5 0.5 0.5 0.25 + 6 0 1 0 0.25 + 7 0 0 1 0.25 + 8 0.5 0.5 0.5 0.25 + + KPOINTS DIRECT_X DIRECT_Y DIRECT_Z WEIGHT + 1 0 0 0 0.25 + 2 0.5 0 0 0.25 + 3 0 0.5 0 0.25 + 4 0.5 0.5 0 0.25 + 5 0 0 0.5 0.25 + 6 0.5 0 0.5 0.25 + 7 0 0.5 0.5 0.25 + 8 0.5 0.5 0.5 0.25 + DONE : INIT K-POINTS Time : 0.0567067414522 (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) = 20 + [fft grid for wave functions] = 24, 24, 24 + number of plane waves = 609 + number of sticks = 91 + + PARALLEL PW FOR WAVE FUNCTIONS + PROC COLUMNS(POT) PW + 1 91 609 + --------------- sum ------------------- + 1 91 609 + DONE : INIT PLANEWAVE Time : 0.0586343016475 (SEC) + + DONE : INIT CHARGE Time : 0.0770398713648 (SEC) + + npwx = 411 + + SETUP NONLOCAL PSEUDOPOTENTIALS IN PLANE WAVE BASIS + Si non-local projectors: + projector 1 L=0 + projector 2 L=0 + projector 3 L=1 + projector 4 L=1 + TOTAL NUMBER OF NONLOCAL PROJECTORS = 16 + DONE : LOCAL POTENTIAL Time : 0.0790216308087 (SEC) + + + Init Non-Local PseudoPotential table : + Init Non-Local-Pseudopotential done. + DONE : NON-LOCAL POTENTIAL Time : 0.0908586103469 (SEC) + + init_chg = atomic + DONE : INIT POTENTIAL Time : 0.119327 (SEC) + + + Make real space PAO into reciprocal space. + max mesh points in Pseudopotential = 601 + dq(describe PAO in reciprocal space) = 0.01 + max q = 542 + + number of pseudo atomic orbitals for Si is 0 + DONE : INIT BASIS Time : 0.188532 (SEC) + + + ------------------------------------------- + STEP OF MOLECULAR DYNAMICS : 0 + ------------------------------------------- + + PW ALGORITHM --------------- ION= 1 ELEC= 1-------------------------------- + + Density error is 0.108050977159 + + PW ALGORITHM --------------- ION= 1 ELEC= 2-------------------------------- + + Density error is 0.00433478793135 + + PW ALGORITHM --------------- ION= 1 ELEC= 3-------------------------------- + + Density error is 0.000134739737693 + + PW ALGORITHM --------------- ION= 1 ELEC= 4-------------------------------- + + Density error is 1.46497067882e-05 + + PW ALGORITHM --------------- ION= 1 ELEC= 5-------------------------------- + + Density error is 5.9940961597e-07 + + PW ALGORITHM --------------- ION= 1 ELEC= 6-------------------------------- + + Density error is 4.85331534904e-08 + + PW ALGORITHM --------------- ION= 1 ELEC= 7-------------------------------- + + Density error is 2.82302502273e-09 + + PW ALGORITHM --------------- ION= 1 ELEC= 8-------------------------------- + + Density error is 6.09165139828e-11 + + charge density convergence is achieved + final etot is -211.771846025 eV + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + Si1 -0.88536267 +0.50046739 +0.15023961 + Si2 +0.88536267 -0.50046739 -0.15023961 + + + ------------------------------------------------------------------------------------------------ + Energy Potential Kinetic Temperature + -7.782469 -7.782469 0 0 + ------------------------------------------------------------------------------------------------ + + + LARGEST GRAD (eV/A) : 0.88536273 + + ------------------------------------------- + STEP OF MOLECULAR DYNAMICS : 1 + ------------------------------------------- + + PW ALGORITHM --------------- ION=2 ELEC=1 -------------------------------- + Notice: Threshold on eigenvalues was too large. + hsover_error=0.08 > DRHO=9.38864e-06 + Origin diag_ethr = 0.01 + New diag_ethr = 1.17358e-07 + + Density error is 2.53667694344e-06 + + PW ALGORITHM --------------- ION=2 ELEC=2 -------------------------------- + + Density error is 2.28596275709e-07 + + PW ALGORITHM --------------- ION=2 ELEC=3 -------------------------------- + + Density error is 1.09611454258e-09 + + PW ALGORITHM --------------- ION=2 ELEC=4 -------------------------------- + + Density error is 1.16027364188e-10 + + charge density convergence is achieved + final etot is -211.781119663 eV + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + Si1 -0.72880473 +0.40857872 +0.10704247 + Si2 +0.72880473 -0.40857872 -0.10704247 + + + ------------------------------------------------------------------------------------------------ + Energy Potential Kinetic Temperature + -7.7824997 -7.7828098 0.0003100856 65.278113 + ------------------------------------------------------------------------------------------------ + + + LARGEST GRAD (eV/A) : 0.72880478 + + ------------------------------------------- + STEP OF MOLECULAR DYNAMICS : 2 + ------------------------------------------- + + PW ALGORITHM --------------- ION=3 ELEC=1 -------------------------------- + Notice: Threshold on eigenvalues was too large. + hsover_error=0.08 > DRHO=5.29133e-05 + Origin diag_ethr = 0.01 + New diag_ethr = 6.61416e-07 + + Density error is 1.76079784755e-05 + + PW ALGORITHM --------------- ION=3 ELEC=2 -------------------------------- + + Density error is 1.54869752683e-06 + + PW ALGORITHM --------------- ION=3 ELEC=3 -------------------------------- + + Density error is 5.21746321277e-09 + + PW ALGORITHM --------------- ION=3 ELEC=4 -------------------------------- + + Density error is 6.06706894933e-10 + + charge density convergence is achieved + final etot is -211.796816626 eV + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + Si1 -0.31603885 +0.17161348 +0.01480266 + Si2 +0.31603885 -0.17161348 -0.01480266 + + + ------------------------------------------------------------------------------------------------ + Energy Potential Kinetic Temperature + -7.7825506 -7.7833867 0.00083609165 176.01103 + ------------------------------------------------------------------------------------------------ + + + LARGEST GRAD (eV/A) : 0.31603887 + + ------------------------------------------- + STEP OF MOLECULAR DYNAMICS : 3 + ------------------------------------------- + + PW ALGORITHM --------------- ION=4 ELEC=1 -------------------------------- + Notice: Threshold on eigenvalues was too large. + hsover_error=0.08 > DRHO=8.44133e-05 + Origin diag_ethr = 0.01 + New diag_ethr = 1.05517e-06 + + Density error is 2.78326201367e-05 + + PW ALGORITHM --------------- ION=4 ELEC=2 -------------------------------- + + Density error is 2.45782918146e-06 + + PW ALGORITHM --------------- ION=4 ELEC=3 -------------------------------- + + Density error is 7.86601506443e-09 + + PW ALGORITHM --------------- ION=4 ELEC=4 -------------------------------- + + Density error is 9.90158000398e-10 + + charge density convergence is achieved + final etot is -211.798755244 eV + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + Si1 +0.20432445 -0.11739811 -0.05295552 + Si2 -0.20432445 +0.11739811 +0.05295552 + + + ------------------------------------------------------------------------------------------------ + Energy Potential Kinetic Temperature + -7.7825562 -7.7834579 0.00090165915 189.81406 + ------------------------------------------------------------------------------------------------ + + + LARGEST GRAD (eV/A) : 0.20432447 + + + -------------------------------------------- + !FINAL_ETOT_IS -211.7987552435292 eV + -------------------------------------------- + + + + + + + |CLASS_NAME---------|NAME---------------|TIME(Sec)-----|CALLS----|AVG------|PER%------- + total 2.7171 17 0.16 1e+02 % + Run_pw plane_wave_line 2.7017 1 2.7 99 % + PW_Basis gathers_scatterp 1.3404 3580 0.00037 49 % + Potential v_of_rho 0.25409 24 0.011 9.4 % + XC_Functional v_xc 0.27163 28 0.0097 10 % + PW_Basis gatherp_scatters 1.3094 3249 0.0004 48 % + Hamilt_PW h_psi 1.7435 5100 0.00034 64 % + Hamilt_PW vloc 1.563 5100 0.00031 58 % + PW_Basis_K recip2real 0.78296 7012 0.00011 29 % + PW_Basis_K real2recip 1.3247 3138 0.00042 49 % + Hamilt_PW vnl 0.14363 5100 2.8e-05 5.3 % + Run_MD_PW md_cells_pw 2.5221 1 2.5 93 % + Run_MD_PW md_ions_pw 2.5221 1 2.5 93 % + FIRE setup 0.94517 1 0.95 35 % + MD_func force_stress 1.5034 2 0.75 55 % + ESolver_KS_PW Run 2.3747 4 0.59 87 % + DiagoCG diag 1.4984 184 0.0081 55 % + DiagoIterAssist diagH_subspace 0.36157 160 0.0023 13 % + MD_func md_force_stress 1.0426 2 0.52 38 % + ---------------------------------------------------------------------------------------- + + CLASS_NAME---------|NAME---------------|MEMORY(MB)-------- + 4.003 + ---------------------------------------------------------- + + Start Time : Thu Jun 9 13:13:06 2022 + Finish Time : Thu Jun 9 13:13:08 2022 + Total Time : 0 h 0 mins 2 secs From f9397124020e9f1dc969137bffc1a1e65c33b937 Mon Sep 17 00:00:00 2001 From: LiuRenxi Date: Sat, 11 Jun 2022 14:02:37 +0800 Subject: [PATCH 4/4] add assertion failed reason for abacus/md interface --- dpdata/abacus/md.py | 8 +-- tests/test_abacus_md.py | 111 +++++++++++++++++++--------------------- 2 files changed, 58 insertions(+), 61 deletions(-) diff --git a/dpdata/abacus/md.py b/dpdata/abacus/md.py index 32774ddca..020f968b5 100644 --- a/dpdata/abacus/md.py +++ b/dpdata/abacus/md.py @@ -37,13 +37,13 @@ def get_coords_from_dump(dumplines, natoms): if "VIRIAL" in dumplines[6]: calc_stress = True else: - assert("POSITIONS" in dumplines[6] and "FORCE" in dumplines[6]) + assert("POSITIONS" in dumplines[6] and "FORCE" in dumplines[6]), "keywords 'POSITIONS' and 'FORCE' cannot be found in the 6th line. Please check." nframes_dump = -1 if calc_stress: nframes_dump = int(nlines/(total_natoms + 13)) else: nframes_dump = int(nlines/(total_natoms + 9)) - assert(nframes_dump > 0) + assert(nframes_dump > 0), "Number of lines in MD_dump file = %d. Number of atoms = %d. The MD_dump file is incomplete."%(nlines, total_natoms) cells = np.zeros([nframes_dump, 3, 3]) stresses = np.zeros([nframes_dump, 3, 3]) forces = np.zeros([nframes_dump, total_natoms, 3]) @@ -66,7 +66,7 @@ def get_coords_from_dump(dumplines, natoms): coords[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+iat])[-6:-3]])*celldm forces[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+iat])[-3:]]) iframe += 1 - assert(iframe == nframes_dump) + assert(iframe == nframes_dump), "iframe=%d, nframe_dump=%d. Number of frames does not match number of lines in MD_dump."%(iframe, nframes_dump) cells *= bohr2ang coords *= bohr2ang stresses *= kbar2evperang3 @@ -80,7 +80,7 @@ def get_energy(outlines, ndump, dump_freq): if nenergy%dump_freq == 0: energy.append(float(line.split()[-2])) nenergy+=1 - assert(ndump == len(energy)) + assert(ndump == len(energy)), "Number of total energies in running_md.log = %d. Number of frames in MD_dump = %d. Please check."%(len(energy), ndump) energy = np.array(energy) return energy diff --git a/tests/test_abacus_md.py b/tests/test_abacus_md.py index 2b40b3840..27fe75c48 100644 --- a/tests/test_abacus_md.py +++ b/tests/test_abacus_md.py @@ -38,67 +38,64 @@ def test_cell(self) : self.assertAlmostEqual(self.system_Si.data['cells'][idx][ii][jj], cell2[ii][jj]) def test_coord(self) : - fp = open('abacus.md/water_coord') - coord = [] - for ii in fp : - coord.append([float(jj) for jj in ii.split()]) - coord = np.array(coord) - coord = coord.reshape([5, 3, 3]) - for ii in range(coord.shape[0]) : - for jj in range(coord.shape[1]) : - for kk in range(coord.shape[2]): - self.assertAlmostEqual(self.system_water.data['coords'][ii][jj][kk], coord[ii][jj][kk]) - fp.close() - - fp2 = open('abacus.md.nostress/Si_coord') - coord = [] - for ii in fp2 : - coord.append([float(jj) for jj in ii.split()]) - coord = np.array(coord) - coord = coord.reshape([4, 2, 3]) - for ii in range(coord.shape[0]) : - for jj in range(coord.shape[1]) : - for kk in range(coord.shape[2]): - self.assertAlmostEqual(self.system_Si.data['coords'][ii][jj][kk], coord[ii][jj][kk]) - fp2.close() + with open('abacus.md/water_coord') as fp: + coord = [] + for ii in fp : + coord.append([float(jj) for jj in ii.split()]) + coord = np.array(coord) + coord = coord.reshape([5, 3, 3]) + for ii in range(coord.shape[0]) : + for jj in range(coord.shape[1]) : + for kk in range(coord.shape[2]): + self.assertAlmostEqual(self.system_water.data['coords'][ii][jj][kk], coord[ii][jj][kk]) + + with open('abacus.md.nostress/Si_coord') as fp2: + coord = [] + for ii in fp2 : + coord.append([float(jj) for jj in ii.split()]) + coord = np.array(coord) + coord = coord.reshape([4, 2, 3]) + for ii in range(coord.shape[0]) : + for jj in range(coord.shape[1]) : + for kk in range(coord.shape[2]): + self.assertAlmostEqual(self.system_Si.data['coords'][ii][jj][kk], coord[ii][jj][kk]) def test_force(self) : - fp = open('abacus.md/water_force') - force = [] - for ii in fp : - force.append([float(jj) for jj in ii.split()]) - force = np.array(force) - force = force.reshape([5, 3, 3]) - for ii in range(force.shape[0]) : - for jj in range(force.shape[1]) : - for kk in range(force.shape[2]): - self.assertAlmostEqual(self.system_water.data['forces'][ii][jj][kk], force[ii][jj][kk]) - fp.close() - - fp2 = open('abacus.md.nostress/Si_force') - force = [] - for ii in fp2 : - force.append([float(jj) for jj in ii.split()]) - force = np.array(force) - force = force.reshape([4, 2, 3]) - for ii in range(force.shape[0]) : - for jj in range(force.shape[1]) : - for kk in range(force.shape[2]): - self.assertAlmostEqual(self.system_Si.data['forces'][ii][jj][kk], force[ii][jj][kk]) - fp2.close() + with open('abacus.md/water_force') as fp: + force = [] + for ii in fp : + force.append([float(jj) for jj in ii.split()]) + force = np.array(force) + force = force.reshape([5, 3, 3]) + for ii in range(force.shape[0]) : + for jj in range(force.shape[1]) : + for kk in range(force.shape[2]): + self.assertAlmostEqual(self.system_water.data['forces'][ii][jj][kk], force[ii][jj][kk]) + + + with open('abacus.md.nostress/Si_force') as fp2: + force = [] + for ii in fp2 : + force.append([float(jj) for jj in ii.split()]) + force = np.array(force) + force = force.reshape([4, 2, 3]) + for ii in range(force.shape[0]) : + for jj in range(force.shape[1]) : + for kk in range(force.shape[2]): + self.assertAlmostEqual(self.system_Si.data['forces'][ii][jj][kk], force[ii][jj][kk]) + def test_virial(self) : - fp = open('abacus.md/water_virial') - virial = [] - for ii in fp : - virial.append([float(jj) for jj in ii.split()]) - virial = np.array(virial) - virial = virial.reshape([5, 3, 3]) - for ii in range(virial.shape[0]) : - for jj in range(virial.shape[1]) : - for kk in range(virial.shape[2]) : - self.assertAlmostEqual(self.system_water.data['virials'][ii][jj][kk], virial[ii][jj][kk]) - fp.close() + with open('abacus.md/water_virial') as fp: + virial = [] + for ii in fp : + virial.append([float(jj) for jj in ii.split()]) + virial = np.array(virial) + virial = virial.reshape([5, 3, 3]) + for ii in range(virial.shape[0]) : + for jj in range(virial.shape[1]) : + for kk in range(virial.shape[2]) : + self.assertAlmostEqual(self.system_water.data['virials'][ii][jj][kk], virial[ii][jj][kk]) def test_energy(self) : ref_energy = np.array([-466.69285117, -466.69929051, -466.69829826, -466.70364664,