Adding a function that could deal with the magnetic_force and spin when using the ABACUS's deltaspin module

.gitattributes

@@ -1 +1,2 @@
 .git_archival.txt  export-subst
+*.py test eol=lf

.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "python.analysis.typeCheckingMode": "basic"
+}

dpdata/abacus/scf.py

@@ -159,6 +159,31 @@
     return Etot, False
 
 
+def get_magnetic(outlines, natoms):
+    Mag = []
+
+    # regular expression
+    mag_pattern = re.compile(
+        r"Total Magnetism on atom:.*\(([\d.\-]+),\s*([\d.\-]+),\s*([\d.\-]+)\)"
+    )
+
+    for line in outlines:
+        # match the magnetic information with regular expression
+        match = mag_pattern.search(line)
+        if match:
+            # Transit to float, and append to the list
+            Mag.append([float(num) for num in match.groups()])
+            # If the lense of the Mag equal to 'natoms', all the magnetic moments collected
+            if len(Mag) == natoms:
+                break
+
+    # If no magnetic information found return an empty NumPy array
+    if len(Mag) == 0:
+        return np.array([[]])
+    else:
+        return np.array(Mag)
+
+
 def collect_force(outlines):
     force = []
     for i, line in enumerate(outlines):
@@ -195,6 +220,42 @@
         return np.array(force[-1])  # only return the last force
 
 
+def collect_mag_force(outlines):
+    mag_force = []
+    for i, line in enumerate(outlines):
+        if "Magnetic force (Ry/uB)" in line:
+            value_pattern = re.compile(
+                r"^\s*ATOM\s+(\d+)\s+[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?\s+[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?\s+[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?\s*$"
+            )
+            j = i
+            # find the first line of force
+            noforce = False
+            while not value_pattern.match(outlines[j]):
+                j += 1
+                if (
+                    j >= i + 10
+                ):  # if can not find the first line of force in 10 lines, then stop
+                    warnings.warn("Warning: can not find the first line of force")
+                    noforce = True
+                    break
+            if noforce:
+                break
+
+            mag_force.append([])
+            while value_pattern.match(outlines[j]):
+                mag_force[-1].append([float(ii) for ii in outlines[j].split()[1:4]])
+                j += 1
+    return mag_force  # only return the last force
+
+
+def get_mag_force(outlines, natoms):
+    mag_force = collect_mag_force(outlines)
+    if len(mag_force) == 0:
+        return [[]]
+    else:
+        return np.array(mag_force[-1])  # only return the last force
+
+
 def collect_stress(outlines):
     stress = []
     for i, line in enumerate(outlines):
@@ -232,6 +293,21 @@
         return np.array(stress[-1]) * kbar2evperang3  # only return the last stress
 
 
+def check_deltaspin(path_in):
+    try:
+        with open(path_in) as file:
+            for line in file:
+                parts = line.strip().split()
+                # Check if the "sc_mag_switch" exist
+                if len(parts) >= 2 and parts[0] == "sc_mag_switch":
+                    return True if parts[1] == "1" else None
+        # if "sc_mag_switch" not exist return None
+        return None
+    except FileNotFoundError:
+        print(f"File not found: {path_in}")
+        return None
+
+
 def get_frame(fname):
     data = {
         "atom_names": [],
@@ -241,6 +317,11 @@
         "coords": [],
         "energies": [],
         "forces": [],
+        "spin": [],
+        "mag_forces": [],
+        "coords_deltaspin": [],
+        "force_deltaspin": [],
+        "deltaspin": [],
     }
 
     if isinstance(fname, str):
@@ -258,6 +339,7 @@
 
     geometry_path_in = get_geometry_in(fname, inlines)
     path_out = get_path_out(fname, inlines)
+    deltaspin = check_deltaspin(path_in)
     if not (CheckFile(geometry_path_in) and CheckFile(path_out)):
         return data
 
@@ -282,10 +364,25 @@
     if stress is not None:
         stress *= np.abs(np.linalg.det(cell))
 
+    if deltaspin is not None:
+        spin = get_magnetic(outlines, natoms)
+        sp_norm = 1.49
+        virtual_len = 0.3
+        mag_forces = get_mag_force(outlines, natoms)
+        coords_deltaspin = np.hstack((coords, coords - spin / sp_norm * virtual_len))
+        force_deltaspin = np.hstack((force, mag_forces))
+
     data["cells"] = cell[np.newaxis, :, :]
     data["coords"] = coords[np.newaxis, :, :]
     data["energies"] = np.array(energy)[np.newaxis]
     data["forces"] = force[np.newaxis, :, :]
+    if deltaspin is not None:
+        data["mag_forces"] = mag_forces[np.newaxis, :, :]
+        data["spin"] = spin[np.newaxis, :, :]
+        data["coords_deltaspin"] = coords_deltaspin[np.newaxis, :, :]
+        data["force_deltaspin"] = force_deltaspin[np.newaxis, :, :]
+        data["deltaspin"] = deltaspin
+
     if stress is not None:
         data["virials"] = stress[np.newaxis, :, :]
     data["orig"] = np.zeros(3)
@@ -297,6 +394,10 @@
     # print("energy = ", data['energies'])
     # print("force = ", data['forces'])
     # print("virial = ", data['virials'])
+    # print("spin = ", data['spin'])
+    # print("mag_forces = ", data['mag_forces'])
+    # print("force_deltaspin = ", data['force_deltaspin'])
+    # print("coords_deltaspin = ", data['coords_deltaspin'])
     return data
 
 

dpdata/deepmd/comp.py

@@ -135,13 +135,25 @@
     nframes = data["cells"].shape[0]
     cells = np.reshape(data["cells"], [nframes, 9]).astype(comp_prec)
     coords = np.reshape(data["coords"], [nframes, -1]).astype(comp_prec)
+    if "coords_deltaspin" in data:
+        coords_deltaspin = np.reshape(data["coords_deltaspin"], [nframes, -1]).astype(
+            comp_prec
+        )
     eners = None
     forces = None
     virials = None
+    force_deltaspin = None
+    deltaspin = None
+    if "deltaspin" in data:
+        deltaspin = data["deltaspin"]
     if "energies" in data:
         eners = np.reshape(data["energies"], [nframes]).astype(comp_prec)
     if "forces" in data:
         forces = np.reshape(data["forces"], [nframes, -1]).astype(comp_prec)
+    if "force_deltaspin" in data:
+        force_deltaspin = np.reshape(data["force_deltaspin"], [nframes, -1]).astype(
+            comp_prec
+        )
     if "virials" in data:
         virials = np.reshape(data["virials"], [nframes, 9]).astype(comp_prec)
     if "atom_pref" in data:
@@ -156,11 +168,21 @@
         set_folder = os.path.join(folder, "set.%03d" % ii)
         os.makedirs(set_folder)
         np.save(os.path.join(set_folder, "box"), cells[set_stt:set_end])
-        np.save(os.path.join(set_folder, "coord"), coords[set_stt:set_end])
+        if deltaspin is not None:
+            np.save(
+                os.path.join(set_folder, "coord"), coords_deltaspin[set_stt:set_end]
+            )
+        else:
+            np.save(os.path.join(set_folder, "coord"), coords[set_stt:set_end])
         if eners is not None:
             np.save(os.path.join(set_folder, "energy"), eners[set_stt:set_end])
         if forces is not None:
-            np.save(os.path.join(set_folder, "force"), forces[set_stt:set_end])
+            if deltaspin is not None:
+                np.save(
+                    os.path.join(set_folder, "force"), force_deltaspin[set_stt:set_end]
+                )
+            else:
+                np.save(os.path.join(set_folder, "force"), forces[set_stt:set_end])
         if virials is not None:
             np.save(os.path.join(set_folder, "virial"), virials[set_stt:set_end])
         if "atom_pref" in data:
@@ -187,6 +209,9 @@
             "energies",
             "forces",
             "virials",
+            "force_deltaspin",
+            "coords_deltaspin",
+            "deltaspin",
         ):
             # skip as these data contains specific rules
             continue

tests/abacus.scf.deltaspin/INPUT

@@ -0,0 +1,29 @@
+INPUT_PARAMETERS
+suffix                          ABACUS
+calculation                     scf
+ecutwfc                         100
+scf_thr                         1.0e-7
+scf_nmax                        200
+out_chg                         0
+smearing_method                 gauss
+smearing_sigma                  0.01
+ks_solver                       genelpa
+basis_type                      lcao
+symmetry                        0
+noncolin			            1
+lspinorb			            0
+nspin                           4
+kspacing                        1.0
+cal_force                       true
+cal_stress                      true
+out_mul                         true
+
+#deltaspin
+sc_mag_switch                   1
+decay_grad_switch               1
+sc_thr                          1e-7
+nsc                             100
+nsc_min                         2
+sc_file                         sc.json
+alpha_trial                     0.01
+sccut                           3