diff --git a/interfaces/ASE_interface/abacuslite/io/generalio.py b/interfaces/ASE_interface/abacuslite/io/generalio.py index 2bf997e6cfe..a2ba99e7f15 100644 --- a/interfaces/ASE_interface/abacuslite/io/generalio.py +++ b/interfaces/ASE_interface/abacuslite/io/generalio.py @@ -220,6 +220,12 @@ def write_stru(stru: Atoms, ind = np.argsort(elem) coords = stru.get_positions()[ind] elem = [elem[i] for i in ind] + + # handle the atomic magnetic moment (issue #6516) + magmoms = np.array([stru[i].magmom for i in ind]).reshape(len(stru), -1) # ncol in [1, 3] + magmoms = [{} if abs(np.linalg.norm(m)) <= 1e-10 + else {'mag': m[0] if len(m) == 1 else ('Cartesian', m.tolist())} + for m in magmoms] elem_uniq, nat = np.unique(elem, return_counts=True) stru_dict = { 'coord_type': 'Cartesian', @@ -236,10 +242,10 @@ def write_stru(stru: Atoms, 'natom': n, 'mag_each': 0.0, 'atom': [ - { - 'coord': coords[j].tolist(), - 'm': [1, 1, 1], - 'v': [0.0, 0.0, 0.0] + magmoms[j] | { + 'coord': coords[j].tolist(), # coordinate + 'm': [1, 1, 1], # mobility + 'v': [0.0, 0.0, 0.0], # velocity } for j in range(np.sum(nat[:i]), np.sum(nat[:i+1])) ] } @@ -745,5 +751,71 @@ def test_load_orbital(self): self.assertEqual(orbital_map['Si'], 'Si_gga_7au_100Ry_2s2p1d.orb') self.assertEqual(orbital_map['In'], 'In_gga_8au_100Ry_2s2p2d1f.orb') + def test_write_stru_with_magmom(self): + # from issue #6516: https://github.com/deepmodeling/abacus-develop/issues/6516 + atoms = Atoms( + symbols=['Co', 'Cr', 'Ni', 'Co'], + cell=np.eye(3), + magmoms=[1, 2, 3, 1] + ) + with tempfile.NamedTemporaryFile(mode='w') as f: + write_stru( + stru=atoms, + outdir=self.testfiles, + pp_file={ # as if we have some pseudopotentials + 'Co': 'Co.upf', 'Cr': 'Cr.upf', 'Ni': 'Ni.upf', 'Co': 'Co.upf', + }, + orb_file={ # as if we have some orbitals + 'Co': 'Co.orb', 'Cr': 'Cr.orb', 'Ni': 'Ni.orb', 'Co': 'Co.orb', + }, + fname=f.name, + ) + stru = read_stru(f.name) + # Co + self.assertEqual(stru['species'][0]['mag_each'], 0.0) + self.assertEqual(stru['species'][0]['atom'][0]['mag'], 1.0) + self.assertEqual(stru['species'][0]['atom'][1]['mag'], 1.0) + # Cr + self.assertEqual(stru['species'][1]['mag_each'], 0.0) + self.assertEqual(stru['species'][1]['atom'][0]['mag'], 2.0) + # Ni + self.assertEqual(stru['species'][2]['mag_each'], 0.0) + self.assertEqual(stru['species'][2]['atom'][0]['mag'], 3.0) + + # then the non-colinear case + atoms = Atoms( + symbols=['Co', 'Cr', 'Ni', 'Co'], + cell=np.eye(3), + magmoms=np.array([ + [1, 0, 0], + [0, 2, 0], + [0, 0, 3], + [1, 0, 0], + ]) + ) + with tempfile.NamedTemporaryFile(mode='w') as f: + write_stru( + stru=atoms, + outdir=self.testfiles, + pp_file={ # as if we have some pseudopotentials + 'Co': 'Co.upf', 'Cr': 'Cr.upf', 'Ni': 'Ni.upf', 'Co': 'Co.upf', + }, + orb_file={ # as if we have some orbitals + 'Co': 'Co.orb', 'Cr': 'Cr.orb', 'Ni': 'Ni.orb', 'Co': 'Co.orb', + }, + fname=f.name, + ) + stru = read_stru(f.name) + # Co + self.assertEqual(stru['species'][0]['mag_each'], 0.0) + self.assertEqual(stru['species'][0]['atom'][0]['mag'], ('Cartesian', [1.0, 0.0, 0.0])) + self.assertEqual(stru['species'][0]['atom'][1]['mag'], ('Cartesian', [1.0, 0.0, 0.0])) + # Cr + self.assertEqual(stru['species'][1]['mag_each'], 0.0) + self.assertEqual(stru['species'][1]['atom'][0]['mag'], ('Cartesian', [0.0, 2.0, 0.0])) + # Ni + self.assertEqual(stru['species'][2]['mag_each'], 0.0) + self.assertEqual(stru['species'][2]['atom'][0]['mag'], ('Cartesian', [0.0, 0.0, 3.0])) + if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/interfaces/ASE_interface/abacuslite/io/latestio.py b/interfaces/ASE_interface/abacuslite/io/latestio.py index 57a38a44299..cba0fce2012 100644 --- a/interfaces/ASE_interface/abacuslite/io/latestio.py +++ b/interfaces/ASE_interface/abacuslite/io/latestio.py @@ -470,7 +470,8 @@ def read_abacus_out(fileobj, ind = ind or list(range(len(frame['elem']))) atoms = Atoms(symbols=np.array(frame['elem'])[ind].tolist(), positions=frame['coords'][ind], - cell=frame['cell']) + cell=frame['cell'], + magmoms=mag[ind]) # from result, a calculator can be assembled # however, sometimes the force and stress is not calculated # in this case, we set them to None diff --git a/interfaces/ASE_interface/abacuslite/io/legacyio.py b/interfaces/ASE_interface/abacuslite/io/legacyio.py index 9772a6363d3..252ff3e3f40 100644 --- a/interfaces/ASE_interface/abacuslite/io/legacyio.py +++ b/interfaces/ASE_interface/abacuslite/io/legacyio.py @@ -662,7 +662,7 @@ def read_magmom_from_running_log(src: str | Path | List[str]) \ while istart < len(raw): ith = None # index of the table header for i, l in enumerate(raw[istart:]): - if re.match(r'\s*Total\sMagnetism\s\(uB\)(\s+x\s+y\s+z)?\s*', l, + if re.match(r'\s*Total\sMagnetism\s\(uB\)(\s+x\s+y\s+z)?\s*$', l, re.IGNORECASE): ith = i break @@ -695,12 +695,12 @@ def read_magmom_from_running_log(src: str | Path | List[str]) \ # the second to fourth items are the x, y, z components res = list(zip(*[l.split() for l in tb_raw]))[1:] assert len(res) in [1, 3] # colinear or non-colinear case - mx, my, mz = (0,) * len(res[-1]), (0,) * len(res[-1]), res[-1] - if len(res) == 3: - mx, my = res[0], res[1] - magmom.append(np.array([list(map(float, (mx, my, mz))) - for mx, my, mz in zip(mx, my, mz)])) - + if len(res) == 1: # colinear case + magmom.append(np.array([list(map(float, res[-1]))]).flatten()) + else: # non-colinear case + mx, my, mz = res + magmom.append(np.array([list(map(float, (mx, my, mz))) + for mx, my, mz in zip(mx, my, mz)])) # update the istart istart += ith + itb + jtb + 1 @@ -828,7 +828,9 @@ def read_abacus_out(fileobj, ind = ind or list(range(len(frame['elem']))) atoms = Atoms(symbols=np.array(frame['elem'])[ind].tolist(), positions=frame['coords'][ind], - cell=frame['cell']) + cell=frame['cell'], + magmoms=mag[ind]) + # from result, a calculator can be assembled # however, sometimes the force and stress is not calculated # in this case, we set them to None @@ -1091,8 +1093,19 @@ def test_read_magmom_from_running_log(self): self.assertEqual(len(magmoms), 2) # 2 cell-relax steps for i, magmom in enumerate(magmoms): self.assertIsInstance(magmom, np.ndarray) - self.assertEqual(magmom.shape, (4, 3)) + self.assertEqual(magmom.shape, (4,)) self.assertAlmostEqual(magmom.sum(), 0.0, delta=1e-3) # AFM + # non-colinear case + fn = self.testfiles / 'pw-symm0-nspin4-gamma-md_' + magmoms = read_magmom_from_running_log(fn) + self.assertIsInstance(magmoms, list) + self.assertEqual(len(magmoms), 2) + for magmom in magmoms: + self.assertTrue( + np.allclose(magmom, + np.array([[0. , 0. , 3.62032142], + [0. , 0. , 3.62032142]]))) + if __name__ == '__main__': unittest.main() diff --git a/interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md b/interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md index fc0508ec7c7..5c30ea1cfdb 100644 --- a/interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md +++ b/interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md @@ -135,6 +135,21 @@ As1 -0.5102760465 -0.5774389693 1.0803355583 Ga1 0.5102760465 0.5774389693 -1.0803355583 ------------------------------------------------------------------------- + + ------------------------------------------------------------------------------------------- + Total Magnetism (uB) x y z +------------------------------------------------------------------------------------------- + As1 0.0000000000 0.0000000000 3.6203214193 + Ga1 0.0000000000 0.0000000000 3.6203214193 +------------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------------- + Total Magnetism (uB) Magnitude (uB) Polar (degree) Azimuth (degree) +------------------------------------------------------------------------------------------- + As1 3.6203214193 0.0000000000 0.0000000000 + Ga1 3.6203214193 0.0000000000 0.0000000000 +------------------------------------------------------------------------------------------- + ------------------------------------------------------------------------------------------------ Energy (Ry) Potential (Ry) Kinetic (Ry) Temperature (K) -142.63058 -142.63343 0.0028501339 300 @@ -194,12 +209,26 @@ As1 0.1191299488 0.2678988387 1.5033416041 Ga1 -0.1191299488 -0.2678988387 -1.5033416041 ------------------------------------------------------------------------- + + ------------------------------------------------------------------------------------------- + Total Magnetism (uB) x y z +------------------------------------------------------------------------------------------- + As1 0.0000000000 0.0000000000 3.6203214193 + Ga1 0.0000000000 0.0000000000 3.6203214193 +------------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------------- + Total Magnetism (uB) Magnitude (uB) Polar (degree) Azimuth (degree) +------------------------------------------------------------------------------------------- + As1 3.6203214193 0.0000000000 0.0000000000 + Ga1 3.6203214193 0.0000000000 0.0000000000 +------------------------------------------------------------------------------------------- + ------------------------------------------------------------------------------------------------ Energy (Ry) Potential (Ry) Kinetic (Ry) Temperature (K) -142.63112 -142.6337 0.0025856408 272.15993 ------------------------------------------------------------------------------------------------ - -------------------------------------------- !FINAL_ETOT_IS -1940.631106372724 eV -------------------------------------------- diff --git a/interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md_ b/interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md_ index 23c2d9a4ae9..45143f8269c 100644 --- a/interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md_ +++ b/interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md_ @@ -159,7 +159,19 @@ total magnetism (Bohr mag/cell) 0 0 1.90445147148 Ga1 -1.0712942948 -2.1781413019 -1.2152201133 ------------------------------------------------------------------------------------------ - +------------------------------------------------------------------------------------------- + Total Magnetism (uB) x y z +------------------------------------------------------------------------------------------- + As1 0.0000000000 0.0000000000 3.6203214193 + Ga1 0.0000000000 0.0000000000 3.6203214193 +------------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------------- + Total Magnetism (uB) Magnitude (uB) Polar (degree) Azimuth (degree) +------------------------------------------------------------------------------------------- + As1 3.6203214193 0.0000000000 0.0000000000 + Ga1 3.6203214193 0.0000000000 0.0000000000 +------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------ Energy (Ry) Potential (Ry) Kinetic (Ry) Temperature (K) @@ -373,16 +385,25 @@ total magnetism (Bohr mag/cell) 0.00000000000 0.00000000000 1.97493083463 Ga1 10.1802174084 5.5311564523 6.8477853698 ------------------------------------------------------------------------------------------ +------------------------------------------------------------------------------------------- + Total Magnetism (uB) x y z +------------------------------------------------------------------------------------------- + As1 0.0000000000 0.0000000000 3.6203214193 + Ga1 0.0000000000 0.0000000000 3.6203214193 +------------------------------------------------------------------------------------------- +------------------------------------------------------------------------------------------- + Total Magnetism (uB) Magnitude (uB) Polar (degree) Azimuth (degree) +------------------------------------------------------------------------------------------- + As1 3.6203214193 0.0000000000 0.0000000000 + Ga1 3.6203214193 0.0000000000 0.0000000000 +------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------ Energy (Ry) Potential (Ry) Kinetic (Ry) Temperature (K) -142.57602 -142.63762 0.061605080 6484.4407 ------------------------------------------------------------------------------------------------ - - - -------------------------------------------- !FINAL_ETOT_IS -1940.684423983896 eV -------------------------------------------- diff --git a/interfaces/ASE_interface/examples/soc.py b/interfaces/ASE_interface/examples/soc.py new file mode 100644 index 00000000000..cdbbd04d6e4 --- /dev/null +++ b/interfaces/ASE_interface/examples/soc.py @@ -0,0 +1,63 @@ +'''this example shows how to perform the noncolinear +spin-orbit coupling calculation''' +import shutil +from pathlib import Path +here = Path(__file__).parent + +pporb = here.parent.parent.parent / 'tests' / 'PP_ORB' + +import numpy as np +from ase.atoms import Atoms +from abacuslite import Abacus, AbacusProfile + +'''SPECIAL: ase does not support the noncolinear spin yet +till 2026/3/24, see ase/outputs.py:L154-155, in which the +magmom cannot be set as the vector, so we release the +datatype of magmom and magmoms by ourself''' +from ase.outputs import _defineprop, all_outputs +del all_outputs['magmom'] +del all_outputs['magmoms'] +_defineprop('magmom', float, shape=3) # re-define the magmom can be set as the vector +_defineprop('magmoms', float, shape=('natoms', 3)) + +fe = Atoms(symbols=['Fe'] * 2, + scaled_positions=[[0., 0., 0.], [0.5, 0.5, 0.5]], + magmoms=np.array([[0, 0, 1], [0, 0, 1]]), + cell=np.eye(3) * 5.2, + pbc=True) + +aprof = AbacusProfile( + command='mpirun -np 2 abacus', + pseudo_dir=pporb, + orbital_dir=pporb, + omp_num_threads=1 +) + +jobdir = here / 'soc' +abacus = Abacus( + profile=aprof, + directory=jobdir, + pseudopotentials={'Fe': 'Fe.upf'}, + basissets={'Fe': 'Fe_gga_6au_100Ry_4s2p2d1f.orb'}, + inp={ + 'basis_type': 'lcao', + 'nspin': 4, + 'lspinorb': True, + 'noncolin': True, + 'out_mul': True, + 'scf_nmax': 10 # not serious setting! + }, + kpts={ + 'mode': 'mp-sampling', + 'gamma-centered': True, + 'nk': (2, 2, 2), + 'kshift': (0, 0, 0) + } +) + +fe.calc = abacus +print(f'Fe (non-colinear): {fe.get_potential_energy()} eV') +for ife, m in enumerate(fe.calc.results['magmoms']): + print(f'Fe{ife}: ({",".join([f"{mi:10.6f}" for mi in m])}) uB') + +shutil.rmtree(jobdir) \ No newline at end of file diff --git a/interfaces/ASE_interface/tests/magnetic.py b/interfaces/ASE_interface/tests/magnetic.py new file mode 100644 index 00000000000..3fbacbdc998 --- /dev/null +++ b/interfaces/ASE_interface/tests/magnetic.py @@ -0,0 +1,129 @@ +''' +test if ABACUS can perform the calculation on the anti-ferromagnetic +and ferromagnetic phases calculation on the BCC Fe +''' +import unittest +import tempfile +import numpy as np +from pathlib import Path +here = Path(__file__).parent +from ase.atoms import Atoms +from abacuslite import AbacusProfile, Abacus + +class TestMagneticSCF(unittest.TestCase): + + def setUp(self): + pporb = here.parent.parent.parent / 'tests' / 'PP_ORB' + self.aprof = AbacusProfile( + command='mpirun -np 2 abacus', + pseudo_dir=pporb, + orbital_dir=pporb, + omp_num_threads=1, + ) + + def test_antiferromagnetic(self): + fe = Atoms(symbols=['Fe'] * 2, + scaled_positions=[[0., 0., 0.], [0.5, 0.5, 0.5]], + magmoms=[-1, 1], + cell=np.eye(3) * 5.2, + pbc=True) + + with tempfile.TemporaryDirectory() as tmpdir: + abacus = Abacus( + profile=self.aprof, + directory=tmpdir, + pseudopotentials={'Fe': 'Fe_ONCV_PBE-1.0.upf'}, + basissets={'Fe': 'Fe_gga_6au_100Ry_4s2p2d1f.orb'}, + inp={ + 'basis_type': 'lcao', + 'gamma_only': True, + 'scf_thr': 1e-3, # fast for test, wrong for production + 'nspin': 2, + 'out_mul': True + } + ) + + fe.calc = abacus + print(f'Fe (AFM): {fe.get_potential_energy()} eV') + for ife, m in enumerate(fe.calc.results['magmoms']): + print(f'Fe{ife}: {m:10.6f} uB') + self.assertGreater(abs(m), 1e-3) + self.assertTrue(np.allclose(np.sum(fe.calc.results['magmoms']), 0.0, atol=1e-6)) + + def test_ferromagnetic(self): + + fe = Atoms(symbols=['Fe'] * 2, + scaled_positions=[[0., 0., 0.], [0.5, 0.5, 0.5]], + magmoms=[1, 1], + cell=np.eye(3) * 5.2, + pbc=True) + + with tempfile.TemporaryDirectory() as tmpdir: + abacus = Abacus( + profile=self.aprof, + directory=tmpdir, + pseudopotentials={'Fe': 'Fe_ONCV_PBE-1.0.upf'}, + basissets={'Fe': 'Fe_gga_6au_100Ry_4s2p2d1f.orb'}, + inp={ + 'basis_type': 'lcao', + 'gamma_only': True, + 'scf_thr': 1e-3, # fast for test, wrong for production + 'nspin': 2, + 'out_mul': True + } + ) + + fe.calc = abacus + print(f'Fe (FM): {fe.get_potential_energy()} eV') + for ife, m in enumerate(fe.calc.results['magmoms']): + print(f'Fe{ife}: {m:10.6f} uB') + self.assertGreater(abs(m), 1e-3) + self.assertGreater(abs(np.sum(fe.calc.results['magmoms'])), 0.0) + + # @unittest.skip('see ase/outputs.py:L154-155, the magmom cannot be set as the vector, so ' + # 'we skip this test for now') + def test_noncolinear(self): + from ase.outputs import _defineprop, all_outputs + # refresh the definition of magmom and magmoms + del all_outputs['magmom'] + del all_outputs['magmoms'] + _defineprop('magmom', float, shape=3) + _defineprop('magmoms', float, shape=('natoms', 3)) + + fe = Atoms(symbols=['Fe'] * 2, + scaled_positions=[[0., 0., 0.], [0.5, 0.5, 0.5]], + magmoms=np.array([[1, 0, 0], [0, 1, 0]]), + cell=np.eye(3) * 5.2, + pbc=True) + + with tempfile.TemporaryDirectory() as tmpdir: + abacus = Abacus( + profile=self.aprof, + directory=tmpdir, + pseudopotentials={'Fe': 'Fe.upf'}, + basissets={'Fe': 'Fe_gga_6au_100Ry_4s2p2d1f.orb'}, + inp={ + 'basis_type': 'lcao', + 'scf_thr': 1e-2, # fast for test, wrong for production + 'nspin': 4, + 'lspinorb': True, + 'noncolin': True, + 'out_mul': True, + 'scf_nmax': 1 + }, + kpts={ + 'mode': 'mp-sampling', + 'gamma-centered': True, + 'nk': (2, 2, 2), + 'kshift': (0, 0, 0) + } + ) + + fe.calc = abacus + print(f'Fe (non-colinear): {fe.get_potential_energy()} eV') + for ife, m in enumerate(fe.calc.results['magmoms']): + print(f'Fe{ife}: ({",".join([f"{mi:10.6f}" for mi in m])}) uB') + self.assertGreater(np.linalg.norm(m), 1e-3) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/interfaces/ASE_interface/tests/xtest.sh b/interfaces/ASE_interface/tests/xtest.sh index 5fbe7fb381b..645c729e7f9 100644 --- a/interfaces/ASE_interface/tests/xtest.sh +++ b/interfaces/ASE_interface/tests/xtest.sh @@ -6,3 +6,4 @@ python3 ./scf.py -v python3 ./relax.py -v python3 ./md.py -v python3 ./band.py -v +python3 ./magnetic.py -v