From cf29c0e031a071de2b1cfd35229027a161409e01 Mon Sep 17 00:00:00 2001 From: root Date: Tue, 5 Nov 2024 14:23:55 +0800 Subject: [PATCH 1/4] abacus: add spin information for fmt=abacus/stru --- dpdata/abacus/md.py | 2 +- dpdata/abacus/relax.py | 2 +- dpdata/abacus/scf.py | 66 +++++++++++++++++++++++++++++++++------ dpdata/plugins/abacus.py | 4 +++ tests/test_abacus_spin.py | 7 +++++ 5 files changed, 70 insertions(+), 11 deletions(-) diff --git a/dpdata/abacus/md.py b/dpdata/abacus/md.py index aa4215d8..f389855f 100644 --- a/dpdata/abacus/md.py +++ b/dpdata/abacus/md.py @@ -167,7 +167,7 @@ def get_frame(fname): with open_file(geometry_path_in) as fp: geometry_inlines = fp.read().split("\n") celldm, cell = get_cell(geometry_inlines) - atom_names, natoms, types, coords, move = get_coords( + atom_names, natoms, types, coords, move, magmom = get_coords( celldm, cell, geometry_inlines, inlines ) # This coords is not to be used. diff --git a/dpdata/abacus/relax.py b/dpdata/abacus/relax.py index 63f13678..ae2fdb13 100644 --- a/dpdata/abacus/relax.py +++ b/dpdata/abacus/relax.py @@ -183,7 +183,7 @@ def get_frame(fname): with open_file(geometry_path_in) as fp: geometry_inlines = fp.read().split("\n") celldm, cell = get_cell(geometry_inlines) - atom_names, natoms, types, coord_tmp, move = get_coords( + atom_names, natoms, types, coord_tmp, move, magmom = get_coords( celldm, cell, geometry_inlines, inlines ) diff --git a/dpdata/abacus/scf.py b/dpdata/abacus/scf.py index 43e65f7e..0406a42a 100644 --- a/dpdata/abacus/scf.py +++ b/dpdata/abacus/scf.py @@ -257,6 +257,49 @@ def parse_stru_pos(pos_line): return pos, move, velocity, magmom, angle1, angle2, constrain, lambda1 +def get_atom_mag_cartesian(atommag,angle1, angle2): + '''Transform atommag, angle1, angle2 to magmom in cartesian coordinates. + + Parameters + ---------- + atommag : float/list of float/None + Atom magnetic moment. + angle1 : float/None + value of angle1. + angle2 : float/None + value of angle2. + + ABACUS support defining mag, angle1, angle2 at the same time. + angle1 is the angle between z-axis and real spin (in degrees). + angle2 is the angle between x-axis and real spin projection in xy-plane (in degrees). + If only mag is defined, then transfer it to magmom directly. + And if mag, angle1, angle2 are defined, then mag is only the norm of magmom, and the direction is defined by angle1 and angle2. + ''' + if atommag is None: + return None + if not (isinstance(atommag, list) or isinstance(atommag, float)): + raise RuntimeError(f"Invalid atommag: {atommag}") + + if angle1 is None and angle2 is None: + if isinstance(atommag, list): + return atommag + else: + return [0,0,atommag] + else: + a1 = 0 + a2 = 0 + if angle1 is not None: + a1 = angle1 + if angle2 is not None: + a2 = angle2 + if isinstance(atommag, list): + mag_norm = np.linalg.norm(atommag) + else: + mag_norm = atommag + return [mag_norm*np.sin(np.radians(a1))*np.cos(np.radians(a2)), + mag_norm*np.sin(np.radians(a1))*np.sin(np.radians(a2)), + mag_norm*np.cos(np.radians(a1))] + def get_coords(celldm, cell, geometry_inlines, inlines=None): coords_lines = get_stru_block(geometry_inlines, "ATOMIC_POSITIONS") @@ -268,9 +311,7 @@ def get_coords(celldm, cell, geometry_inlines, inlines=None): coords = [] # coordinations of atoms move = [] # move flag of each atom velocity = [] # velocity of each atom - mag = [] # magnetic moment of each atom - angle1 = [] # angle1 of each atom - angle2 = [] # angle2 of each atom + mags = [] # magnetic moment of each atom sc = [] # spin constraint flag of each atom lambda_ = [] # lambda of each atom @@ -278,6 +319,7 @@ def get_coords(celldm, cell, geometry_inlines, inlines=None): line_idx = 1 # starting line of first element for it in range(ntype): atom_names.append(coords_lines[line_idx].split()[0]) + atom_type_mag = float(coords_lines[line_idx+1].split()[0]) line_idx += 2 atom_numbs.append(int(coords_lines[line_idx].split()[0])) line_idx += 1 @@ -302,17 +344,20 @@ def get_coords(celldm, cell, geometry_inlines, inlines=None): if imove is not None: move.append(imove) velocity.append(ivelocity) - mag.append(imagmom) - angle1.append(iangle1) - angle2.append(iangle2) sc.append(iconstrain) lambda_.append(ilambda1) + + # calculate the magnetic moment in cartesian coordinates + mag = get_atom_mag_cartesian(imagmom, iangle1, iangle2) + if mag is None: + mag = [0,0,atom_type_mag] + mags.append(mag) line_idx += 1 coords = np.array(coords) # need transformation!!! atom_types = np.array(atom_types) move = np.array(move, dtype=bool) - return atom_names, atom_numbs, atom_types, coords, move + return atom_names, atom_numbs, atom_types, coords, move, mags def get_energy(outlines): @@ -479,9 +524,11 @@ def get_frame(fname): outlines = fp.read().split("\n") celldm, cell = get_cell(geometry_inlines) - atom_names, natoms, types, coords, move = get_coords( + atom_names, natoms, types, coords, move, magmom = get_coords( # here the magmom is the initial magnetic moment in STRU celldm, cell, geometry_inlines, inlines ) + + magmom, magforce = get_mag_force(outlines) if len(magmom) > 0: magmom = magmom[-1:] @@ -565,7 +612,7 @@ def get_frame_from_stru(fname): nele = get_nele_from_stru(geometry_inlines) inlines = [f"ntype {nele}"] celldm, cell = get_cell(geometry_inlines) - atom_names, natoms, types, coords, move = get_coords( + atom_names, natoms, types, coords, move, magmom = get_coords( celldm, cell, geometry_inlines, inlines ) data = {} @@ -575,6 +622,7 @@ def get_frame_from_stru(fname): data["cells"] = cell[np.newaxis, :, :] data["coords"] = coords[np.newaxis, :, :] data["orig"] = np.zeros(3) + data["spins"] = np.array([magmom]) if len(move) > 0: data["move"] = move[np.newaxis, :, :] diff --git a/dpdata/plugins/abacus.py b/dpdata/plugins/abacus.py index b3e7c98a..f97ab846 100644 --- a/dpdata/plugins/abacus.py +++ b/dpdata/plugins/abacus.py @@ -20,6 +20,8 @@ @Format.register("stru") class AbacusSTRUFormat(Format): def from_system(self, file_name, **kwargs): + data = dpdata.abacus.scf.get_frame(file_name) + register_mag_data(data) return dpdata.abacus.scf.get_frame_from_stru(file_name) def to_system(self, data, file_name: FileType, frame_idx=0, **kwargs): @@ -55,6 +57,7 @@ def register_mag_data(data): required=False, deepmd_name="spin", ) + dpdata.System.register_data_type(dt) dpdata.LabeledSystem.register_data_type(dt) if "mag_forces" in data: dt = DataType( @@ -64,6 +67,7 @@ def register_mag_data(data): required=False, deepmd_name="force_mag", ) + dpdata.System.register_data_type(dt) dpdata.LabeledSystem.register_data_type(dt) diff --git a/tests/test_abacus_spin.py b/tests/test_abacus_spin.py index df8ef7d0..a786eb7b 100644 --- a/tests/test_abacus_spin.py +++ b/tests/test_abacus_spin.py @@ -155,3 +155,10 @@ def test_md(self): np.testing.assert_almost_equal( data["mag_forces"], sys2.data["mag_forces"], decimal=8 ) + + def test_read_stru_spin(self): + mysys = dpdata.System("abacus.spin/STRU", fmt="abacus/stru") + self.assertTrue("spins" in mysys.data) + np.testing.assert_almost_equal(mysys.data["spins"],[[[0,0,2.4],[0,0,2.4]]],decimal=8) + + From 46592e9c811e424e53690e67375e9481f220da48 Mon Sep 17 00:00:00 2001 From: root Date: Tue, 5 Nov 2024 15:53:55 +0800 Subject: [PATCH 2/4] add test --- tests/abacus.spin/STRU.spin | 24 ++++++++++++++++++++++++ tests/test_abacus_spin.py | 17 +++++++++++++++-- 2 files changed, 39 insertions(+), 2 deletions(-) create mode 100644 tests/abacus.spin/STRU.spin diff --git a/tests/abacus.spin/STRU.spin b/tests/abacus.spin/STRU.spin new file mode 100644 index 00000000..340430c9 --- /dev/null +++ b/tests/abacus.spin/STRU.spin @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +Fe 55.845 Fe.upf + +NUMERICAL_ORBITAL +Fe_gga_10au_200.0Ry_4s2p2d1f.orb + +LATTICE_CONSTANT +1.880277359 +LATTICE_VECTORS + 2.8274254848 0.0000000000 0.0000000000 #latvec1 + 0.0000000000 2.8274254848 0.0000000000 #latvec2 + 0.0000000000 0.0000000000 2.8274254848 #latvec3 + +ATOMIC_POSITIONS +Direct + +Fe #label +1 #magnetism +4 #number of atoms + 0.0000000000 0.000000000 0.000000000 mag 0 0 2 + 0.1000000000 0.1000000000 0.1000000000 mag 3 + 0.2000000000 0.2000000000 0.2000000000 mag 3 angle1 90 + 0.3000000000 0.3000000000 0.3000000000 mag 3 4 0 angle1 90 angle2 90 + diff --git a/tests/test_abacus_spin.py b/tests/test_abacus_spin.py index a786eb7b..65761d77 100644 --- a/tests/test_abacus_spin.py +++ b/tests/test_abacus_spin.py @@ -157,8 +157,21 @@ def test_md(self): ) def test_read_stru_spin(self): - mysys = dpdata.System("abacus.spin/STRU", fmt="abacus/stru") + mysys = dpdata.System("abacus.spin/STRU.spin", fmt="abacus/stru") self.assertTrue("spins" in mysys.data) - np.testing.assert_almost_equal(mysys.data["spins"],[[[0,0,2.4],[0,0,2.4]]],decimal=8) + print(mysys.data["spins"]) + + """ + 0.0000000000 0.000000000 0.000000000 mag 0 0 2 + 0.1000000000 0.1000000000 0.1000000000 mag 3 + 0.2000000000 0.2000000000 0.2000000000 mag 3 angle1 90 + 0.3000000000 0.3000000000 0.3000000000 mag 3 4 0 angle1 90 angle2 90 + """ + np.testing.assert_almost_equal(mysys.data["spins"][0][0],[0,0,2],decimal=8) + np.testing.assert_almost_equal(mysys.data["spins"][0][1],[0,0,3],decimal=8) + np.testing.assert_almost_equal(mysys.data["spins"][0][2],[3,0,0],decimal=8) + np.testing.assert_almost_equal(mysys.data["spins"][0][3],[0,5,0],decimal=8) + + From 222148fbbfd0bb8bdc9366827883addde0417bd8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 5 Nov 2024 07:56:04 +0000 Subject: [PATCH 3/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- dpdata/abacus/scf.py | 43 +++++++++++++++++++++------------------ tests/test_abacus_spin.py | 16 ++++++--------- 2 files changed, 29 insertions(+), 30 deletions(-) diff --git a/dpdata/abacus/scf.py b/dpdata/abacus/scf.py index 0406a42a..d489e702 100644 --- a/dpdata/abacus/scf.py +++ b/dpdata/abacus/scf.py @@ -257,9 +257,10 @@ def parse_stru_pos(pos_line): return pos, move, velocity, magmom, angle1, angle2, constrain, lambda1 -def get_atom_mag_cartesian(atommag,angle1, angle2): - '''Transform atommag, angle1, angle2 to magmom in cartesian coordinates. - + +def get_atom_mag_cartesian(atommag, angle1, angle2): + """Transform atommag, angle1, angle2 to magmom in cartesian coordinates. + Parameters ---------- atommag : float/list of float/None @@ -267,24 +268,23 @@ def get_atom_mag_cartesian(atommag,angle1, angle2): angle1 : float/None value of angle1. angle2 : float/None - value of angle2. - + value of angle2. ABACUS support defining mag, angle1, angle2 at the same time. angle1 is the angle between z-axis and real spin (in degrees). - angle2 is the angle between x-axis and real spin projection in xy-plane (in degrees). + angle2 is the angle between x-axis and real spin projection in xy-plane (in degrees). If only mag is defined, then transfer it to magmom directly. And if mag, angle1, angle2 are defined, then mag is only the norm of magmom, and the direction is defined by angle1 and angle2. - ''' + """ if atommag is None: return None if not (isinstance(atommag, list) or isinstance(atommag, float)): raise RuntimeError(f"Invalid atommag: {atommag}") - + if angle1 is None and angle2 is None: if isinstance(atommag, list): return atommag else: - return [0,0,atommag] + return [0, 0, atommag] else: a1 = 0 a2 = 0 @@ -296,10 +296,12 @@ def get_atom_mag_cartesian(atommag,angle1, angle2): mag_norm = np.linalg.norm(atommag) else: mag_norm = atommag - return [mag_norm*np.sin(np.radians(a1))*np.cos(np.radians(a2)), - mag_norm*np.sin(np.radians(a1))*np.sin(np.radians(a2)), - mag_norm*np.cos(np.radians(a1))] - + return [ + mag_norm * np.sin(np.radians(a1)) * np.cos(np.radians(a2)), + mag_norm * np.sin(np.radians(a1)) * np.sin(np.radians(a2)), + mag_norm * np.cos(np.radians(a1)), + ] + def get_coords(celldm, cell, geometry_inlines, inlines=None): coords_lines = get_stru_block(geometry_inlines, "ATOMIC_POSITIONS") @@ -319,7 +321,7 @@ def get_coords(celldm, cell, geometry_inlines, inlines=None): line_idx = 1 # starting line of first element for it in range(ntype): atom_names.append(coords_lines[line_idx].split()[0]) - atom_type_mag = float(coords_lines[line_idx+1].split()[0]) + atom_type_mag = float(coords_lines[line_idx + 1].split()[0]) line_idx += 2 atom_numbs.append(int(coords_lines[line_idx].split()[0])) line_idx += 1 @@ -346,11 +348,11 @@ def get_coords(celldm, cell, geometry_inlines, inlines=None): velocity.append(ivelocity) sc.append(iconstrain) lambda_.append(ilambda1) - + # calculate the magnetic moment in cartesian coordinates mag = get_atom_mag_cartesian(imagmom, iangle1, iangle2) if mag is None: - mag = [0,0,atom_type_mag] + mag = [0, 0, atom_type_mag] mags.append(mag) line_idx += 1 @@ -524,11 +526,12 @@ def get_frame(fname): outlines = fp.read().split("\n") celldm, cell = get_cell(geometry_inlines) - atom_names, natoms, types, coords, move, magmom = get_coords( # here the magmom is the initial magnetic moment in STRU - celldm, cell, geometry_inlines, inlines + atom_names, natoms, types, coords, move, magmom = ( + get_coords( # here the magmom is the initial magnetic moment in STRU + celldm, cell, geometry_inlines, inlines + ) ) - - + magmom, magforce = get_mag_force(outlines) if len(magmom) > 0: magmom = magmom[-1:] diff --git a/tests/test_abacus_spin.py b/tests/test_abacus_spin.py index 65761d77..7f2bec83 100644 --- a/tests/test_abacus_spin.py +++ b/tests/test_abacus_spin.py @@ -155,23 +155,19 @@ def test_md(self): np.testing.assert_almost_equal( data["mag_forces"], sys2.data["mag_forces"], decimal=8 ) - + def test_read_stru_spin(self): mysys = dpdata.System("abacus.spin/STRU.spin", fmt="abacus/stru") self.assertTrue("spins" in mysys.data) print(mysys.data["spins"]) - + """ 0.0000000000 0.000000000 0.000000000 mag 0 0 2 0.1000000000 0.1000000000 0.1000000000 mag 3 0.2000000000 0.2000000000 0.2000000000 mag 3 angle1 90 0.3000000000 0.3000000000 0.3000000000 mag 3 4 0 angle1 90 angle2 90 """ - np.testing.assert_almost_equal(mysys.data["spins"][0][0],[0,0,2],decimal=8) - np.testing.assert_almost_equal(mysys.data["spins"][0][1],[0,0,3],decimal=8) - np.testing.assert_almost_equal(mysys.data["spins"][0][2],[3,0,0],decimal=8) - np.testing.assert_almost_equal(mysys.data["spins"][0][3],[0,5,0],decimal=8) - - - - + np.testing.assert_almost_equal(mysys.data["spins"][0][0], [0, 0, 2], decimal=8) + np.testing.assert_almost_equal(mysys.data["spins"][0][1], [0, 0, 3], decimal=8) + np.testing.assert_almost_equal(mysys.data["spins"][0][2], [3, 0, 0], decimal=8) + np.testing.assert_almost_equal(mysys.data["spins"][0][3], [0, 5, 0], decimal=8) From 35e0d833a3fe8db350a4de95a7a2fdca03d35921 Mon Sep 17 00:00:00 2001 From: root Date: Wed, 6 Nov 2024 14:52:46 +0800 Subject: [PATCH 4/4] fix --- dpdata/plugins/abacus.py | 4 ++-- tests/abacus.scf/stru_test | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/dpdata/plugins/abacus.py b/dpdata/plugins/abacus.py index 79ba3a4e..b40d9966 100644 --- a/dpdata/plugins/abacus.py +++ b/dpdata/plugins/abacus.py @@ -20,9 +20,9 @@ @Format.register("stru") class AbacusSTRUFormat(Format): def from_system(self, file_name, **kwargs): - data = dpdata.abacus.scf.get_frame(file_name) + data = dpdata.abacus.scf.get_frame_from_stru(file_name) register_mag_data(data) - return dpdata.abacus.scf.get_frame_from_stru(file_name) + return data def to_system(self, data, file_name: FileType, frame_idx=0, **kwargs): """Dump the system into ABACUS STRU format file. diff --git a/tests/abacus.scf/stru_test b/tests/abacus.scf/stru_test index e4036409..c3a1917d 100644 --- a/tests/abacus.scf/stru_test +++ b/tests/abacus.scf/stru_test @@ -22,11 +22,11 @@ Cartesian # Cartesian(Unit is LATTICE_CONSTANT) C 0.0 1 -5.192682633809 4.557725978258 4.436846615358 1 1 1 +5.192682633809 4.557725978258 4.436846615358 1 1 1 mag 0.000000000000 0.000000000000 0.000000000000 H 0.0 4 -5.416431453540 4.011298860305 3.511161492417 0 0 0 -4.131588222365 4.706745191323 4.431136645083 1 0 1 -5.630930319126 5.521640894956 4.450356541303 1 0 1 -5.499851012568 4.003388899277 5.342621842622 0 1 1 +5.416431453540 4.011298860305 3.511161492417 0 0 0 mag 0.000000000000 0.000000000000 0.000000000000 +4.131588222365 4.706745191323 4.431136645083 1 0 1 mag 0.000000000000 0.000000000000 0.000000000000 +5.630930319126 5.521640894956 4.450356541303 1 0 1 mag 0.000000000000 0.000000000000 0.000000000000 +5.499851012568 4.003388899277 5.342621842622 0 1 1 mag 0.000000000000 0.000000000000 0.000000000000