Skip to content

Commit e97cf59

Browse files
kirk0830kirk0830
authored
Feature: support the setting on atomic magnetic moment via ase.atoms.Atoms.magmom (#7124)
* Feature: support the setting on atomic magnetic moment via ase.atoms.Atoms.magmom * add tests * update the testfiles * add test for non-colinear case
1 parent 86e9314 commit e97cf59

8 files changed

Lines changed: 348 additions & 19 deletions

File tree

‎interfaces/ASE_interface/abacuslite/io/generalio.py‎

Lines changed: 76 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -220,6 +220,12 @@ def write_stru(stru: Atoms,
220220
ind = np.argsort(elem)
221221
coords = stru.get_positions()[ind]
222222
elem = [elem[i] for i in ind]
223+
224+
# handle the atomic magnetic moment (issue #6516)
225+
magmoms = np.array([stru[i].magmom for i in ind]).reshape(len(stru), -1) # ncol in [1, 3]
226+
magmoms = [{} if abs(np.linalg.norm(m)) <= 1e-10
227+
else {'mag': m[0] if len(m) == 1 else ('Cartesian', m.tolist())}
228+
for m in magmoms]
223229
elem_uniq, nat = np.unique(elem, return_counts=True)
224230
stru_dict = {
225231
'coord_type': 'Cartesian',
@@ -236,10 +242,10 @@ def write_stru(stru: Atoms,
236242
'natom': n,
237243
'mag_each': 0.0,
238244
'atom': [
239-
{
240-
'coord': coords[j].tolist(),
241-
'm': [1, 1, 1],
242-
'v': [0.0, 0.0, 0.0]
245+
magmoms[j] | {
246+
'coord': coords[j].tolist(), # coordinate
247+
'm': [1, 1, 1], # mobility
248+
'v': [0.0, 0.0, 0.0], # velocity
243249
} for j in range(np.sum(nat[:i]), np.sum(nat[:i+1]))
244250
]
245251
}
@@ -745,5 +751,71 @@ def test_load_orbital(self):
745751
self.assertEqual(orbital_map['Si'], 'Si_gga_7au_100Ry_2s2p1d.orb')
746752
self.assertEqual(orbital_map['In'], 'In_gga_8au_100Ry_2s2p2d1f.orb')
747753

754+
def test_write_stru_with_magmom(self):
755+
# from issue #6516: https://github.com/deepmodeling/abacus-develop/issues/6516
756+
atoms = Atoms(
757+
symbols=['Co', 'Cr', 'Ni', 'Co'],
758+
cell=np.eye(3),
759+
magmoms=[1, 2, 3, 1]
760+
)
761+
with tempfile.NamedTemporaryFile(mode='w') as f:
762+
write_stru(
763+
stru=atoms,
764+
outdir=self.testfiles,
765+
pp_file={ # as if we have some pseudopotentials
766+
'Co': 'Co.upf', 'Cr': 'Cr.upf', 'Ni': 'Ni.upf', 'Co': 'Co.upf',
767+
},
768+
orb_file={ # as if we have some orbitals
769+
'Co': 'Co.orb', 'Cr': 'Cr.orb', 'Ni': 'Ni.orb', 'Co': 'Co.orb',
770+
},
771+
fname=f.name,
772+
)
773+
stru = read_stru(f.name)
774+
# Co
775+
self.assertEqual(stru['species'][0]['mag_each'], 0.0)
776+
self.assertEqual(stru['species'][0]['atom'][0]['mag'], 1.0)
777+
self.assertEqual(stru['species'][0]['atom'][1]['mag'], 1.0)
778+
# Cr
779+
self.assertEqual(stru['species'][1]['mag_each'], 0.0)
780+
self.assertEqual(stru['species'][1]['atom'][0]['mag'], 2.0)
781+
# Ni
782+
self.assertEqual(stru['species'][2]['mag_each'], 0.0)
783+
self.assertEqual(stru['species'][2]['atom'][0]['mag'], 3.0)
784+
785+
# then the non-colinear case
786+
atoms = Atoms(
787+
symbols=['Co', 'Cr', 'Ni', 'Co'],
788+
cell=np.eye(3),
789+
magmoms=np.array([
790+
[1, 0, 0],
791+
[0, 2, 0],
792+
[0, 0, 3],
793+
[1, 0, 0],
794+
])
795+
)
796+
with tempfile.NamedTemporaryFile(mode='w') as f:
797+
write_stru(
798+
stru=atoms,
799+
outdir=self.testfiles,
800+
pp_file={ # as if we have some pseudopotentials
801+
'Co': 'Co.upf', 'Cr': 'Cr.upf', 'Ni': 'Ni.upf', 'Co': 'Co.upf',
802+
},
803+
orb_file={ # as if we have some orbitals
804+
'Co': 'Co.orb', 'Cr': 'Cr.orb', 'Ni': 'Ni.orb', 'Co': 'Co.orb',
805+
},
806+
fname=f.name,
807+
)
808+
stru = read_stru(f.name)
809+
# Co
810+
self.assertEqual(stru['species'][0]['mag_each'], 0.0)
811+
self.assertEqual(stru['species'][0]['atom'][0]['mag'], ('Cartesian', [1.0, 0.0, 0.0]))
812+
self.assertEqual(stru['species'][0]['atom'][1]['mag'], ('Cartesian', [1.0, 0.0, 0.0]))
813+
# Cr
814+
self.assertEqual(stru['species'][1]['mag_each'], 0.0)
815+
self.assertEqual(stru['species'][1]['atom'][0]['mag'], ('Cartesian', [0.0, 2.0, 0.0]))
816+
# Ni
817+
self.assertEqual(stru['species'][2]['mag_each'], 0.0)
818+
self.assertEqual(stru['species'][2]['atom'][0]['mag'], ('Cartesian', [0.0, 0.0, 3.0]))
819+
748820
if __name__ == '__main__':
749821
unittest.main()

‎interfaces/ASE_interface/abacuslite/io/latestio.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -470,7 +470,8 @@ def read_abacus_out(fileobj,
470470
ind = ind or list(range(len(frame['elem'])))
471471
atoms = Atoms(symbols=np.array(frame['elem'])[ind].tolist(),
472472
positions=frame['coords'][ind],
473-
cell=frame['cell'])
473+
cell=frame['cell'],
474+
magmoms=mag[ind])
474475
# from result, a calculator can be assembled
475476
# however, sometimes the force and stress is not calculated
476477
# in this case, we set them to None

‎interfaces/ASE_interface/abacuslite/io/legacyio.py‎

Lines changed: 22 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -662,7 +662,7 @@ def read_magmom_from_running_log(src: str | Path | List[str]) \
662662
while istart < len(raw):
663663
ith = None # index of the table header
664664
for i, l in enumerate(raw[istart:]):
665-
if re.match(r'\s*Total\sMagnetism\s\(uB\)(\s+x\s+y\s+z)?\s*', l,
665+
if re.match(r'\s*Total\sMagnetism\s\(uB\)(\s+x\s+y\s+z)?\s*$', l,
666666
re.IGNORECASE):
667667
ith = i
668668
break
@@ -695,12 +695,12 @@ def read_magmom_from_running_log(src: str | Path | List[str]) \
695695
# the second to fourth items are the x, y, z components
696696
res = list(zip(*[l.split() for l in tb_raw]))[1:]
697697
assert len(res) in [1, 3] # colinear or non-colinear case
698-
mx, my, mz = (0,) * len(res[-1]), (0,) * len(res[-1]), res[-1]
699-
if len(res) == 3:
700-
mx, my = res[0], res[1]
701-
magmom.append(np.array([list(map(float, (mx, my, mz)))
702-
for mx, my, mz in zip(mx, my, mz)]))
703-
698+
if len(res) == 1: # colinear case
699+
magmom.append(np.array([list(map(float, res[-1]))]).flatten())
700+
else: # non-colinear case
701+
mx, my, mz = res
702+
magmom.append(np.array([list(map(float, (mx, my, mz)))
703+
for mx, my, mz in zip(mx, my, mz)]))
704704
# update the istart
705705
istart += ith + itb + jtb + 1
706706

@@ -828,7 +828,9 @@ def read_abacus_out(fileobj,
828828
ind = ind or list(range(len(frame['elem'])))
829829
atoms = Atoms(symbols=np.array(frame['elem'])[ind].tolist(),
830830
positions=frame['coords'][ind],
831-
cell=frame['cell'])
831+
cell=frame['cell'],
832+
magmoms=mag[ind])
833+
832834
# from result, a calculator can be assembled
833835
# however, sometimes the force and stress is not calculated
834836
# in this case, we set them to None
@@ -1091,8 +1093,19 @@ def test_read_magmom_from_running_log(self):
10911093
self.assertEqual(len(magmoms), 2) # 2 cell-relax steps
10921094
for i, magmom in enumerate(magmoms):
10931095
self.assertIsInstance(magmom, np.ndarray)
1094-
self.assertEqual(magmom.shape, (4, 3))
1096+
self.assertEqual(magmom.shape, (4,))
10951097
self.assertAlmostEqual(magmom.sum(), 0.0, delta=1e-3) # AFM
10961098

1099+
# non-colinear case
1100+
fn = self.testfiles / 'pw-symm0-nspin4-gamma-md_'
1101+
magmoms = read_magmom_from_running_log(fn)
1102+
self.assertIsInstance(magmoms, list)
1103+
self.assertEqual(len(magmoms), 2)
1104+
for magmom in magmoms:
1105+
self.assertTrue(
1106+
np.allclose(magmom,
1107+
np.array([[0. , 0. , 3.62032142],
1108+
[0. , 0. , 3.62032142]])))
1109+
10971110
if __name__ == '__main__':
10981111
unittest.main()

‎interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md‎

Lines changed: 30 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -135,6 +135,21 @@
135135
As1 -0.5102760465 -0.5774389693 1.0803355583
136136
Ga1 0.5102760465 0.5774389693 -1.0803355583
137137
-------------------------------------------------------------------------
138+
139+
-------------------------------------------------------------------------------------------
140+
Total Magnetism (uB) x y z
141+
-------------------------------------------------------------------------------------------
142+
As1 0.0000000000 0.0000000000 3.6203214193
143+
Ga1 0.0000000000 0.0000000000 3.6203214193
144+
-------------------------------------------------------------------------------------------
145+
146+
-------------------------------------------------------------------------------------------
147+
Total Magnetism (uB) Magnitude (uB) Polar (degree) Azimuth (degree)
148+
-------------------------------------------------------------------------------------------
149+
As1 3.6203214193 0.0000000000 0.0000000000
150+
Ga1 3.6203214193 0.0000000000 0.0000000000
151+
-------------------------------------------------------------------------------------------
152+
138153
------------------------------------------------------------------------------------------------
139154
Energy (Ry) Potential (Ry) Kinetic (Ry) Temperature (K)
140155
-142.63058 -142.63343 0.0028501339 300
@@ -194,12 +209,26 @@
194209
As1 0.1191299488 0.2678988387 1.5033416041
195210
Ga1 -0.1191299488 -0.2678988387 -1.5033416041
196211
-------------------------------------------------------------------------
212+
213+
-------------------------------------------------------------------------------------------
214+
Total Magnetism (uB) x y z
215+
-------------------------------------------------------------------------------------------
216+
As1 0.0000000000 0.0000000000 3.6203214193
217+
Ga1 0.0000000000 0.0000000000 3.6203214193
218+
-------------------------------------------------------------------------------------------
219+
220+
-------------------------------------------------------------------------------------------
221+
Total Magnetism (uB) Magnitude (uB) Polar (degree) Azimuth (degree)
222+
-------------------------------------------------------------------------------------------
223+
As1 3.6203214193 0.0000000000 0.0000000000
224+
Ga1 3.6203214193 0.0000000000 0.0000000000
225+
-------------------------------------------------------------------------------------------
226+
197227
------------------------------------------------------------------------------------------------
198228
Energy (Ry) Potential (Ry) Kinetic (Ry) Temperature (K)
199229
-142.63112 -142.6337 0.0025856408 272.15993
200230
------------------------------------------------------------------------------------------------
201231

202-
203232
--------------------------------------------
204233
!FINAL_ETOT_IS -1940.631106372724 eV
205234
--------------------------------------------

‎interfaces/ASE_interface/abacuslite/io/testfiles/pw-symm0-nspin4-gamma-md_‎

Lines changed: 25 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -159,7 +159,19 @@ total magnetism (Bohr mag/cell) 0 0 1.90445147148
159159
Ga1 -1.0712942948 -2.1781413019 -1.2152201133
160160
------------------------------------------------------------------------------------------
161161

162-
162+
-------------------------------------------------------------------------------------------
163+
Total Magnetism (uB) x y z
164+
-------------------------------------------------------------------------------------------
165+
As1 0.0000000000 0.0000000000 3.6203214193
166+
Ga1 0.0000000000 0.0000000000 3.6203214193
167+
-------------------------------------------------------------------------------------------
168+
169+
-------------------------------------------------------------------------------------------
170+
Total Magnetism (uB) Magnitude (uB) Polar (degree) Azimuth (degree)
171+
-------------------------------------------------------------------------------------------
172+
As1 3.6203214193 0.0000000000 0.0000000000
173+
Ga1 3.6203214193 0.0000000000 0.0000000000
174+
-------------------------------------------------------------------------------------------
163175

164176
------------------------------------------------------------------------------------------------
165177
Energy (Ry) Potential (Ry) Kinetic (Ry) Temperature (K)
@@ -373,16 +385,25 @@ total magnetism (Bohr mag/cell) 0.00000000000 0.00000000000 1.97493083463
373385
Ga1 10.1802174084 5.5311564523 6.8477853698
374386
------------------------------------------------------------------------------------------
375387

388+
-------------------------------------------------------------------------------------------
389+
Total Magnetism (uB) x y z
390+
-------------------------------------------------------------------------------------------
391+
As1 0.0000000000 0.0000000000 3.6203214193
392+
Ga1 0.0000000000 0.0000000000 3.6203214193
393+
-------------------------------------------------------------------------------------------
376394

395+
-------------------------------------------------------------------------------------------
396+
Total Magnetism (uB) Magnitude (uB) Polar (degree) Azimuth (degree)
397+
-------------------------------------------------------------------------------------------
398+
As1 3.6203214193 0.0000000000 0.0000000000
399+
Ga1 3.6203214193 0.0000000000 0.0000000000
400+
-------------------------------------------------------------------------------------------
377401

378402
------------------------------------------------------------------------------------------------
379403
Energy (Ry) Potential (Ry) Kinetic (Ry) Temperature (K)
380404
-142.57602 -142.63762 0.061605080 6484.4407
381405
------------------------------------------------------------------------------------------------
382406

383-
384-
385-
386407
--------------------------------------------
387408
!FINAL_ETOT_IS -1940.684423983896 eV
388409
--------------------------------------------

‎interfaces/ASE_interface/examples/soc.py‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
'''this example shows how to perform the noncolinear
2+
spin-orbit coupling calculation'''
3+
import shutil
4+
from pathlib import Path
5+
here = Path(__file__).parent
6+
7+
pporb = here.parent.parent.parent / 'tests' / 'PP_ORB'
8+
9+
import numpy as np
10+
from ase.atoms import Atoms
11+
from abacuslite import Abacus, AbacusProfile
12+
13+
'''SPECIAL: ase does not support the noncolinear spin yet
14+
till 2026/3/24, see ase/outputs.py:L154-155, in which the
15+
magmom cannot be set as the vector, so we release the
16+
datatype of magmom and magmoms by ourself'''
17+
from ase.outputs import _defineprop, all_outputs
18+
del all_outputs['magmom']
19+
del all_outputs['magmoms']
20+
_defineprop('magmom', float, shape=3) # re-define the magmom can be set as the vector
21+
_defineprop('magmoms', float, shape=('natoms', 3))
22+
23+
fe = Atoms(symbols=['Fe'] * 2,
24+
scaled_positions=[[0., 0., 0.], [0.5, 0.5, 0.5]],
25+
magmoms=np.array([[0, 0, 1], [0, 0, 1]]),
26+
cell=np.eye(3) * 5.2,
27+
pbc=True)
28+
29+
aprof = AbacusProfile(
30+
command='mpirun -np 2 abacus',
31+
pseudo_dir=pporb,
32+
orbital_dir=pporb,
33+
omp_num_threads=1
34+
)
35+
36+
jobdir = here / 'soc'
37+
abacus = Abacus(
38+
profile=aprof,
39+
directory=jobdir,
40+
pseudopotentials={'Fe': 'Fe.upf'},
41+
basissets={'Fe': 'Fe_gga_6au_100Ry_4s2p2d1f.orb'},
42+
inp={
43+
'basis_type': 'lcao',
44+
'nspin': 4,
45+
'lspinorb': True,
46+
'noncolin': True,
47+
'out_mul': True,
48+
'scf_nmax': 10 # not serious setting!
49+
},
50+
kpts={
51+
'mode': 'mp-sampling',
52+
'gamma-centered': True,
53+
'nk': (2, 2, 2),
54+
'kshift': (0, 0, 0)
55+
}
56+
)
57+
58+
fe.calc = abacus
59+
print(f'Fe (non-colinear): {fe.get_potential_energy()} eV')
60+
for ife, m in enumerate(fe.calc.results['magmoms']):
61+
print(f'Fe{ife}: ({",".join([f"{mi:10.6f}" for mi in m])}) uB')
62+
63+
shutil.rmtree(jobdir)

0 commit comments

Comments
 (0)