diff --git a/dpdata/deepmd/comp.py b/dpdata/deepmd/comp.py index 67eddd9f..b7a259e0 100644 --- a/dpdata/deepmd/comp.py +++ b/dpdata/deepmd/comp.py @@ -26,7 +26,8 @@ def _load_set(folder, nopbc: bool): cells = np.zeros((coords.shape[0], 3, 3)) else: cells = np.load(os.path.join(folder, "box.npy")) - return cells, coords + spins = _cond_load_data(os.path.join(folder, "spin.npy")) + return cells, coords, spins def to_system_data(folder, type_map=None, labels=True): @@ -38,13 +39,18 @@ def to_system_data(folder, type_map=None, labels=True): sets = sorted(glob.glob(os.path.join(folder, "set.*"))) all_cells = [] all_coords = [] + all_spins = [] for ii in sets: - cells, coords = _load_set(ii, data.get("nopbc", False)) + cells, coords, spins = _load_set(ii, data.get("nopbc", False)) nframes = np.reshape(cells, [-1, 3, 3]).shape[0] all_cells.append(np.reshape(cells, [nframes, 3, 3])) all_coords.append(np.reshape(coords, [nframes, -1, 3])) + if spins is not None: + all_spins.append(np.reshape(spins, [nframes, -1, 3])) data["cells"] = np.concatenate(all_cells, axis=0) data["coords"] = np.concatenate(all_coords, axis=0) + if len(all_spins) > 0: + data["spins"] = np.concatenate(all_spins, axis=0) # allow custom dtypes if labels: dtypes = dpdata.system.LabeledSystem.DTYPES @@ -59,6 +65,7 @@ def to_system_data(folder, type_map=None, labels=True): "orig", "cells", "coords", + "spins", "real_atom_names", "nopbc", ): diff --git a/dpdata/deepmd/raw.py b/dpdata/deepmd/raw.py index 81b01cc0..e7cdb1f1 100644 --- a/dpdata/deepmd/raw.py +++ b/dpdata/deepmd/raw.py @@ -64,6 +64,7 @@ def to_system_data(folder, type_map=None, labels=True): "orig", "cells", "coords", + "spins", "real_atom_types", "real_atom_names", "nopbc", diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py index 72fee27d..d3d0e4a8 100644 --- a/dpdata/lammps/dump.py +++ b/dpdata/lammps/dump.py @@ -132,6 +132,42 @@ def safe_get_posi(lines, cell, orig=np.zeros(3), unwrap=False): ) # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions +def get_spintype(keys): + key_sp = ["sp", "spx", "spy", "spz"] + key_csp = ["c_spin[1]", "c_spin[2]", "c_spin[3]", "c_spin[4]"] + lmp_sp_type = [key_sp, key_csp] + for k in range(2): + if all(i in keys for i in lmp_sp_type[k]): + return lmp_sp_type[k] + + +def safe_get_spin_force(lines): + blk, head = _get_block(lines, "ATOMS") + keys = head.split() + sp_type = get_spintype(keys) + assert sp_type is not None, "Dump file does not contain spin!" + id_idx = keys.index("id") - 2 + sp = keys.index(sp_type[0]) - 2 + spx = keys.index(sp_type[1]) - 2 + spy = keys.index(sp_type[2]) - 2 + spz = keys.index(sp_type[3]) - 2 + sp_force = [] + for ii in blk: + words = ii.split() + sp_force.append( + [ + float(words[id_idx]), + float(words[sp]), + float(words[spx]), + float(words[spy]), + float(words[spz]), + ] + ) + sp_force.sort() + sp_force = np.array(sp_force)[:, 1:] + return sp_force + + def get_dumpbox(lines): blk, h = _get_block(lines, "BOX BOUNDS") bounds = np.zeros([3, 2]) @@ -216,6 +252,12 @@ def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False): system["cells"] = [np.array(cell)] system["atom_types"] = get_atype(lines, type_idx_zero=type_idx_zero) system["coords"] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] + contain_spin = False + blk, head = _get_block(lines, "ATOMS") + if "sp" in head: + contain_spin = True + spin_force = safe_get_spin_force(lines) + system["spins"] = [spin_force[:, :1] * spin_force[:, 1:4]] for ii in range(1, len(array_lines)): bounds, tilt = get_dumpbox(array_lines[ii]) orig, cell = dumpbox2box(bounds, tilt) @@ -228,6 +270,13 @@ def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False): system["coords"].append( safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap)[idx] ) + if contain_spin: + system["spins"].append( + safe_get_spin_force(array_lines[ii])[:, :1] + * safe_get_spin_force(array_lines[ii])[:, 1:4] + ) + if contain_spin: + system["spins"] = np.array(system["spins"]) system["cells"] = np.array(system["cells"]) system["coords"] = np.array(system["coords"]) return system diff --git a/dpdata/lammps/lmp.py b/dpdata/lammps/lmp.py index 604b18d1..e78443c3 100644 --- a/dpdata/lammps/lmp.py +++ b/dpdata/lammps/lmp.py @@ -127,6 +127,19 @@ def get_posi(lines): return np.array(posis) +def get_spins(lines): + atom_lines = get_atoms(lines) + if len(atom_lines[0].split()) < 8: + return None + spins_ori = [] + spins_norm = [] + for ii in atom_lines: + spins_ori.append([float(jj) for jj in ii.split()[5:8]]) + spins_norm.append([float(jj) for jj in ii.split()[-1:]]) + spins = np.array(spins_ori) * np.array(spins_norm) + return spins + + def get_lmpbox(lines): box_info = [] tilt = np.zeros(3) @@ -151,6 +164,9 @@ def get_lmpbox(lines): def system_data(lines, type_map=None, type_idx_zero=True): system = {} system["atom_numbs"] = get_natoms_vec(lines) + spins = get_spins(lines) + if spins is not None: + system["spins"] = np.array([spins]) system["atom_names"] = [] if type_map is None: for ii in range(len(system["atom_numbs"])): @@ -216,14 +232,54 @@ def from_system_data(system, f_idx=0): + ptr_float_fmt + "\n" ) - for ii in range(natoms): - ret += coord_fmt % ( - ii + 1, - system["atom_types"][ii] + 1, - system["coords"][f_idx][ii][0] - system["orig"][0], - system["coords"][f_idx][ii][1] - system["orig"][1], - system["coords"][f_idx][ii][2] - system["orig"][2], + if "spins" in system.keys(): + coord_fmt = ( + coord_fmt.strip("\n") + + " " + + ptr_float_fmt + + " " + + ptr_float_fmt + + " " + + ptr_float_fmt + + " " + + ptr_float_fmt + + "\n" ) + spins_norm = np.linalg.norm(system["spins"][f_idx], axis=1) + for ii in range(natoms): + if "spins" in system.keys(): + if spins_norm[ii] != 0: + ret += coord_fmt % ( + ii + 1, + system["atom_types"][ii] + 1, + system["coords"][f_idx][ii][0] - system["orig"][0], + system["coords"][f_idx][ii][1] - system["orig"][1], + system["coords"][f_idx][ii][2] - system["orig"][2], + system["spins"][f_idx][ii][0] / spins_norm[ii], + system["spins"][f_idx][ii][1] / spins_norm[ii], + system["spins"][f_idx][ii][2] / spins_norm[ii], + spins_norm[ii], + ) + else: + ret += coord_fmt % ( + ii + 1, + system["atom_types"][ii] + 1, + system["coords"][f_idx][ii][0] - system["orig"][0], + system["coords"][f_idx][ii][1] - system["orig"][1], + system["coords"][f_idx][ii][2] - system["orig"][2], + system["spins"][f_idx][ii][0], + system["spins"][f_idx][ii][1], + system["spins"][f_idx][ii][2] + 1, + spins_norm[ii], + ) + else: + ret += coord_fmt % ( + ii + 1, + system["atom_types"][ii] + 1, + system["coords"][f_idx][ii][0] - system["orig"][0], + system["coords"][f_idx][ii][1] - system["orig"][1], + system["coords"][f_idx][ii][2] - system["orig"][2], + ) return ret diff --git a/dpdata/plugins/vasp_deltaspin.py b/dpdata/plugins/vasp_deltaspin.py new file mode 100644 index 00000000..2a9d74df --- /dev/null +++ b/dpdata/plugins/vasp_deltaspin.py @@ -0,0 +1,107 @@ +from __future__ import annotations + +import re + +import numpy as np + +import dpdata.vasp_deltaspin.outcar +import dpdata.vasp_deltaspin.poscar +from dpdata.format import Format +from dpdata.utils import uniq_atom_names + + +@Format.register("vasp_deltaspin/poscar") +@Format.register("vasp_deltaspin/contcar") +class VASPPoscarFormat(Format): + @Format.post("rot_lower_triangular") + def from_system(self, file_name, **kwargs): + with open(file_name) as fp: + lines = [line.rstrip("\n") for line in fp] + with open(file_name[:-6] + "INCAR") as fp: + lines_incar = [line.rstrip("\n") for line in fp] + data = dpdata.vasp_deltaspin.poscar.to_system_data(lines, lines_incar) + data = uniq_atom_names(data) + return data + + def to_system(self, data, file_name, frame_idx=0, **kwargs): + """Dump the system in vasp POSCAR format. + + Parameters + ---------- + data : dict + The system data + file_name : str + The output file name + frame_idx : int + The index of the frame to dump + **kwargs : dict + other parameters + """ + w_str, m_str = VASPStringFormat().to_system(data, frame_idx=frame_idx) + with open(file_name, "w") as fp: + fp.write(w_str) + + with open(file_name[:-6] + "INCAR") as fp: + tmp_incar = fp.read() + res_incar = re.sub( + r"MAGMOM[\s\S]*?\n\nM_CONST[\s\S]*?\n\n", m_str, tmp_incar, re.S + ) + with open(file_name[:-6] + "INCAR", "w") as fp: + fp.write(res_incar) + + +@Format.register("vasp/string") +class VASPStringFormat(Format): + def to_system(self, data, frame_idx=0, **kwargs): + """Dump the system in vasp POSCAR format string. + + Parameters + ---------- + data : dict + The system data + frame_idx : int + The index of the frame to dump + **kwargs : dict + other parameters + """ + assert frame_idx < len(data["coords"]) + return dpdata.vasp_deltaspin.poscar.from_system_data(data, frame_idx) + + +# rotate the system to lammps convention +@Format.register("vasp_deltaspin/outcar") +class VASPOutcarFormat(Format): + @Format.post("rot_lower_triangular") + def from_labeled_system( + self, file_name, begin=0, step=1, convergence_check=True, **kwargs + ): + data = {} + ml = kwargs.get("ml", False) + ( + data["atom_names"], + data["atom_numbs"], + data["atom_types"], + data["cells"], + data["coords"], + data["spins"], + data["energies"], + data["forces"], + data["mag_forces"], + tmp_virial, + ) = dpdata.vasp_deltaspin.outcar.get_frames( + file_name, + begin=begin, + step=step, + ml=ml, + convergence_check=convergence_check, + ) + if tmp_virial is not None: + data["virials"] = tmp_virial + # scale virial to the unit of eV + if "virials" in data: + v_pref = 1 * 1e3 / 1.602176621e6 + for ii in range(data["cells"].shape[0]): + vol = np.linalg.det(np.reshape(data["cells"][ii], [3, 3])) + data["virials"][ii] *= v_pref * vol + data = uniq_atom_names(data) + return data diff --git a/dpdata/system.py b/dpdata/system.py index 08136cb9..2f4414f3 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -95,6 +95,7 @@ class System: DataType( "coords", np.ndarray, (Axis.NFRAMES, Axis.NATOMS, 3), deepmd_name="coord" ), + DataType("spins", np.ndarray, (Axis.NFRAMES, Axis.NATOMS, 3), required=False), DataType( "real_atom_types", np.ndarray, (Axis.NFRAMES, Axis.NATOMS), required=False ), @@ -712,6 +713,10 @@ def affine_map(self, trans, f_idx: int | numbers.Integral = 0): assert np.linalg.det(trans) != 0 self.data["cells"][f_idx] = np.matmul(self.data["cells"][f_idx], trans) self.data["coords"][f_idx] = np.matmul(self.data["coords"][f_idx], trans) + try: + self.data["spins"][f_idx] = np.matmul(self.data["spins"][f_idx], trans) + except: + pass @post_funcs.register("shift_orig_zero") def _shift_orig_zero(self): @@ -1210,6 +1215,9 @@ class LabeledSystem(System): DataType( "forces", np.ndarray, (Axis.NFRAMES, Axis.NATOMS, 3), deepmd_name="force" ), + DataType( + "mag_forces", np.ndarray, (Axis.NFRAMES, Axis.NATOMS, 3), required=False + ), DataType( "virials", np.ndarray, @@ -1793,3 +1801,5 @@ def to_format(self, *args, **kwargs): add_format_methods() + +# %% diff --git a/dpdata/vasp_deltaspin/__init__.py b/dpdata/vasp_deltaspin/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/dpdata/vasp_deltaspin/outcar.py b/dpdata/vasp_deltaspin/outcar.py new file mode 100644 index 00000000..192106b1 --- /dev/null +++ b/dpdata/vasp_deltaspin/outcar.py @@ -0,0 +1,199 @@ +from __future__ import annotations + +import re +import warnings + +import numpy as np + + +def system_info(lines, type_idx_zero=False): + atom_names = [] + atom_numbs = None + nelm = None + for ii in lines: + ii_word_list = ii.split() + if "TITEL" in ii: + # get atom names from POTCAR info, tested only for PAW_PBE ... + _ii = ii.split()[3] + if "_" in _ii: + # for case like : TITEL = PAW_PBE Sn_d 06Sep2000 + atom_names.append(_ii.split("_")[0]) + else: + atom_names.append(_ii) + # a stricker check for "NELM"; compatible with distingct formats in different versions(6 and older, newers_expect-to-work) of vasp + elif nelm is None: + m = re.search(r"NELM\s*=\s*(\d+)", ii) + if m: + nelm = int(m.group(1)) + if "ions per type" in ii: + atom_numbs_ = [int(s) for s in ii.split()[4:]] + if atom_numbs is None: + atom_numbs = atom_numbs_ + else: + assert atom_numbs == atom_numbs_, "in consistent numb atoms in OUTCAR" + assert nelm is not None, "cannot find maximum steps for each SC iteration" + assert atom_numbs is not None, "cannot find ion type info in OUTCAR" + atom_names = atom_names[: len(atom_numbs)] + atom_types = [] + for idx, ii in enumerate(atom_numbs): + for jj in range(ii): + if type_idx_zero: + atom_types.append(idx) + else: + atom_types.append(idx + 1) + return atom_names, atom_numbs, np.array(atom_types, dtype=int), nelm + + +def get_outcar_block(fp, ml=False): + blk = [] + energy_token = ["free energy TOTEN", "free energy ML TOTEN"] + ml_index = int(ml) + for ii in fp: + if not ii: + return blk + blk.append(ii.rstrip("\n")) + if energy_token[ml_index] in ii: + return blk + return blk + + +# we assume that the force is printed ... +def get_frames(fname, begin=0, step=1, ml=False, convergence_check=True): + fp = open(fname) + blk = get_outcar_block(fp) + + atom_names, atom_numbs, atom_types, nelm = system_info(blk, type_idx_zero=True) + ntot = sum(atom_numbs) + + all_coords = [] + all_spins = [] + all_cells = [] + all_energies = [] + all_forces = [] + all_mag_forces = [] + all_virials = [] + + cc = 0 + rec_failed = [] + while len(blk) > 0: + if cc >= begin and (cc - begin) % step == 0: + coord, cell, energy, force, virial, is_converge = analyze_block( + blk, ntot, nelm, ml + ) + if len(coord) == 0: + break + if is_converge or not convergence_check: + all_coords.append(coord) + all_cells.append(cell) + all_energies.append(energy) + all_forces.append(force) + if virial is not None: + all_virials.append(virial) + if not is_converge: + rec_failed.append(cc + 1) + + blk = get_outcar_block(fp, ml) + cc += 1 + + with open(fname[:-6] + "OSZICAR") as f: + oszicar_blk = f.read() + mag_blk = re.findall(r"Magnetic Force(.*?)\n \n", oszicar_blk, re.S)[-1].split( + "\n" + ) + mag_force = [list(map(float, ss.split()[-3:])) for ss in mag_blk[1:]] + spin_blk = re.findall(r"MW_current(.*?)\n ion", oszicar_blk, re.S)[-1].split("\n") + spin = [list(map(float, ss.split()[1:4])) for ss in spin_blk[1:]] + all_mag_forces.append(mag_force) + all_spins.append(spin) + + if len(rec_failed) > 0: + prt = ( + "so they are not collected." + if convergence_check + else "but they are still collected due to the requirement for ignoring convergence checks." + ) + warnings.warn( + f"The following structures were unconverged: {rec_failed}; " + prt + ) + + if len(all_virials) == 0: + all_virials = None + else: + all_virials = np.array(all_virials) + fp.close() + return ( + atom_names, + atom_numbs, + atom_types, + np.array(all_cells), + np.array(all_coords), + np.array(all_spins), + np.array(all_energies), + np.array(all_forces), + np.array(all_mag_forces), + all_virials, + ) + + +def analyze_block(lines, ntot, nelm, ml=False): + coord = [] + cell = [] + energy = None + force = [] + virial = None + is_converge = True + sc_index = 0 + # select different searching tokens based on the ml label + energy_token = ["free energy TOTEN", "free energy ML TOTEN"] + energy_index = [4, 5] + virial_token = ["FORCE on cell =-STRESS in cart. coord. units", "ML FORCE"] + virial_index = [14, 4] + cell_token = ["VOLUME and BASIS", "ML FORCE"] + cell_index = [5, 12] + ml_index = int(ml) + for idx, ii in enumerate(lines): + # if set ml == True, is_converged will always be True + if ("Iteration" in ii) and (not ml): + sc_index = int(ii.split()[3][:-1]) + if sc_index >= nelm: + is_converge = False + elif energy_token[ml_index] in ii: + energy = float(ii.split()[energy_index[ml_index]]) + if len(force) == 0: + raise ValueError("cannot find forces in OUTCAR block") + if len(coord) == 0: + raise ValueError("cannot find coordinates in OUTCAR block") + if len(cell) == 0: + raise ValueError("cannot find cell in OUTCAR block") + return coord, cell, energy, force, virial, is_converge + elif cell_token[ml_index] in ii: + for dd in range(3): + tmp_l = lines[idx + cell_index[ml_index] + dd] + cell.append([float(ss) for ss in tmp_l.replace("-", " -").split()[0:3]]) + elif virial_token[ml_index] in ii: + in_kB_index = virial_index[ml_index] + while idx + in_kB_index < len(lines) and ( + not lines[idx + in_kB_index].split()[0:2] == ["in", "kB"] + ): + in_kB_index += 1 + assert idx + in_kB_index < len( + lines + ), 'ERROR: "in kB" is not found in OUTCAR. Unable to extract virial.' + tmp_v = [float(ss) for ss in lines[idx + in_kB_index].split()[2:8]] + virial = np.zeros([3, 3]) + virial[0][0] = tmp_v[0] + virial[1][1] = tmp_v[1] + virial[2][2] = tmp_v[2] + virial[0][1] = tmp_v[3] + virial[1][0] = tmp_v[3] + virial[1][2] = tmp_v[4] + virial[2][1] = tmp_v[4] + virial[0][2] = tmp_v[5] + virial[2][0] = tmp_v[5] + elif "TOTAL-FORCE" in ii and (("ML" in ii) == ml): + for jj in range(idx + 2, idx + 2 + ntot): + tmp_l = lines[jj] + info = [float(ss) for ss in tmp_l.split()] + coord.append(info[:3]) + force.append(info[3:6]) + return coord, cell, energy, force, virial, is_converge diff --git a/dpdata/vasp_deltaspin/poscar.py b/dpdata/vasp_deltaspin/poscar.py new file mode 100644 index 00000000..7c7905fb --- /dev/null +++ b/dpdata/vasp_deltaspin/poscar.py @@ -0,0 +1,120 @@ +#!/usr/bin/python3 +from __future__ import annotations + +import numpy as np + + +def _to_system_data_lower(lines, lines_incar, cartesian=True): + """Treat as cartesian poscar.""" + system = {} + system["atom_names"] = [str(ii) for ii in lines[5].split()] + system["atom_numbs"] = [int(ii) for ii in lines[6].split()] + scale = float(lines[1]) + cell = [] + for ii in range(2, 5): + boxv = [float(jj) for jj in lines[ii].split()] + boxv = np.array(boxv) * scale + cell.append(boxv) + system["cells"] = [np.array(cell)] + natoms = sum(system["atom_numbs"]) + coord = [] + for ii in range(8, 8 + natoms): + tmpv = [float(jj) for jj in lines[ii].split()[:3]] + if cartesian: + tmpv = np.array(tmpv) * scale + else: + tmpv = np.matmul(np.array(tmpv), system["cells"][0]) + coord.append(tmpv) + system["coords"] = [np.array(coord)] + system["orig"] = np.zeros(3) + atom_types = [] + for idx, ii in enumerate(system["atom_numbs"]): + for jj in range(ii): + atom_types.append(idx) + system["atom_types"] = np.array(atom_types, dtype=int) + system["cells"] = np.array(system["cells"]) + system["coords"] = np.array(system["coords"]) + + for idx, incar_param in enumerate(lines_incar): + if "MAGMOM" in incar_param: + start_index = idx + elif "M_CONSTR" in incar_param: + end_index = idx + system["spins"] = [lines_incar[start_index].replace("=", "").strip().split()[1:4]] + for idx in range(start_index + 1, end_index - 1): + system["spins"].append(lines_incar[idx].strip().split()[:3]) + system["spins"] = np.array([system["spins"]]).astype("float64") + count = np.sum(np.linalg.norm(system["spins"][0], axis=1) > 0) + spins_idx = np.where(np.cumsum(system["atom_numbs"]) <= count) + return system + + +def to_system_data(lines, lines_incar): + # remove the line that has 'selective dynamics' + if lines[7][0] == "S" or lines[7][0] == "s": + lines.pop(7) + is_cartesian = lines[7][0] in ["C", "c", "K", "k"] + if not is_cartesian: + if lines[7][0] not in ["d", "D"]: + raise RuntimeError( + "seem not to be a valid POSCAR of vasp 5.x, may be a POSCAR of vasp 4.x?" + ) + return _to_system_data_lower(lines, lines_incar, is_cartesian) + + +def from_system_data(system, f_idx=0, skip_zeros=True): + ret = "" + for ii, name in zip(system["atom_numbs"], system["atom_names"]): + if ii == 0: + continue + ret += "%s%d " % (name, ii) + ret += "\n" + ret += "1.0\n" + for ii in system["cells"][f_idx]: + for jj in ii: + ret += "%.16e " % jj + ret += "\n" + for idx, ii in enumerate(system["atom_names"]): + if system["atom_numbs"][idx] == 0: + continue + ret += "%s " % ii + ret += "\n" + for ii in system["atom_numbs"]: + if ii == 0: + continue + ret += "%d " % ii + ret += "\n" + # should use Cartesian for VESTA software + ret += "Cartesian\n" + atype = system["atom_types"] + posis = system["coords"][f_idx] + # atype_idx = [[idx,tt] for idx,tt in enumerate(atype)] + # sort_idx = np.argsort(atype, kind = 'mergesort') + sort_idx = np.lexsort((np.arange(len(atype)), atype)) + atype = atype[sort_idx] + posis = posis[sort_idx] + posi_list = [] + for ii in posis: + posi_list.append(f"{ii[0]:15.10f} {ii[1]:15.10f} {ii[2]:15.10f}") + posi_list.append("") + ret += "\n".join(posi_list) + + magmom_incar = "MAGMOM = " + mconstr_incar = "M_CONSTR = " + for idx, tmp_spin in enumerate(system["spins"][f_idx]): + if idx == 0: + magmom_incar += ( + f"{tmp_spin[0]:10.5f} {tmp_spin[1]:10.5f} {tmp_spin[2]:10.5f} \ \n" + ) + mconstr_incar += ( + f"{tmp_spin[0]:10.5f} {tmp_spin[1]:10.5f} {tmp_spin[2]:10.5f} \ \n" + ) + elif idx == len(system["spins"][f_idx]) - 1: + magmom_incar += f" {tmp_spin[0]:10.5f} {tmp_spin[1]:10.5f} {tmp_spin[2]:10.5f} \n" + mconstr_incar += f" {tmp_spin[0]:10.5f} {tmp_spin[1]:10.5f} {tmp_spin[2]:10.5f} \n" + else: + magmom_incar += f" {tmp_spin[0]:10.5f} {tmp_spin[1]:10.5f} {tmp_spin[2]:10.5f} \ \n" + mconstr_incar += f" {tmp_spin[0]:10.5f} {tmp_spin[1]:10.5f} {tmp_spin[2]:10.5f} \ \n" + magmom_incar += "\n" + mconstr_incar + "\n" + + return ret, magmom_incar