Improve robustness of openmx/md cell parsing and support new output format

dpdata/openmx/omx.py

@@ -27,23 +27,6 @@
 import warnings
 from collections import OrderedDict
 
-### iterout.c from OpenMX soure code: column numbers and physical quantities ###
-# /* 1: */
-# /* 2,3,4: */
-# /* 5,6,7: force *
-# /* 8: x-component of velocity */
-# /* 9: y-component of velocity */
-# /* 10: z-component of velocity */
-# /* 11: Net charge, electron charge is defined to be negative. */
-# /* 12: magnetic moment (muB) */
-# /* 13,14: angles of spin */
-
-# 15: scf_convergence_flag (optional)
-#
-# 1. Move the declaration of `scf_convergence_flag` in `DFT.c` to `openmx_common.h`.
-# 2. Add `scf_convergence_flag` output to the end of `iterout.c` where `*.md` is written.
-# 3. Recompile OpenMX.
-
 
 def load_atom(lines):
     atom_names = []
@@ -56,9 +39,8 @@ def load_atom(lines):
         elif atom_names_mode:
             parts = line.split()
             atom_names.append(parts[1])
-    natoms = len(atom_names)
     atom_names_original = atom_names
-    atom_names = list(OrderedDict.fromkeys(set(atom_names)))  # Python>=3.7
+    atom_names = list(OrderedDict.fromkeys(set(atom_names)))
     atom_names = sorted(
         atom_names, key=atom_names_original.index
     )  # Unique ordering of atomic species
@@ -82,24 +64,20 @@ def load_atom(lines):
 
 
 def load_cells(lines):
-    cell, cells = [], []
-    for index, line in enumerate(lines):
+    cells = []
+    for line in lines:
         if "Cell_Vectors=" in line:
-            parts = line.split()
-            if len(parts) == 21:  # MD.Type is NVT_NH
-                cell.append([float(parts[12]), float(parts[13]), float(parts[14])])
-                cell.append([float(parts[15]), float(parts[16]), float(parts[17])])
-                cell.append([float(parts[18]), float(parts[19]), float(parts[20])])
-            elif len(parts) == 16:  # MD.Type is Opt
-                cell.append([float(parts[7]), float(parts[8]), float(parts[9])])
-                cell.append([float(parts[10]), float(parts[11]), float(parts[12])])
-                cell.append([float(parts[13]), float(parts[14]), float(parts[15])])
-            else:
-                raise RuntimeError(
-                    "Does the file System.Name.md contain unsupported calculation results?"
-                )
+            part = line.split("Cell_Vectors=")[1]
+            parts = part.split()
+            values = list(map(float, parts[:9]))
+            cell = [values[0:3], values[3:6], values[6:9]]
             cells.append(cell)
-            cell = []
+            # Checking SCF converged or not
+            for token in line.split():
+                if token.startswith("scf_conv="):
+                    scf_conv = int(token.split("=")[1])
+                    if scf_conv == 0:
+                        warnings.warn("SCF not converged!", stacklevel=2)
     cells = np.array(cells)
     return cells
 
@@ -119,7 +97,7 @@ def load_param_file(fname: FileType, mdname: FileType):
 def load_coords(lines, atom_names, natoms):
     cnt = 0
     coord, coords = [], []
-    for index, line in enumerate(lines):
+    for line in lines:
         if "time=" in line:
             continue
         for atom_name in atom_names:
@@ -129,9 +107,6 @@ def load_coords(lines, atom_names, natoms):
                 parts = line.split()
                 for_line = [float(parts[1]), float(parts[2]), float(parts[3])]
                 coord.append(for_line)
-                # It may be necessary to recompile OpenMX to make scf convergence determination.
-                if len(parts) == 15 and parts[14] == "0":
-                    warnings.warn("SCF in System.Name.md has not converged!")
         if cnt == natoms:
             coords.append(coord)
             cnt = 0
@@ -180,7 +155,7 @@ def load_energy(lines):
 def load_force(lines, atom_names, atom_numbs):
     cnt = 0
     field, fields = [], []
-    for index, line in enumerate(lines):
+    for line in lines:
         if "time=" in line:
             continue
         for atom_name in atom_names:
@@ -209,7 +184,7 @@ def to_system_label(fname, mdname):
 
 
 if __name__ == "__main__":
-    file_name = "Cdia"
+    file_name = "Au111Surface"
     fname = f"{file_name}.dat"
     mdname = f"{file_name}.md"
     atom_names, atom_numbs, atom_types, cells = load_param_file(fname, mdname)
@@ -219,7 +194,7 @@ def to_system_label(fname, mdname):
     print(atom_names)
     print(atom_numbs)
     print(atom_types)
-    # print(cells.shape)
-    # print(coords.shape)
-    # print(len(energy))
-    # print(force.shape)
+# print(cells.shape)
+# print(coords.shape)
+# print(len(energy))
+# print(force.shape)

tests/openmx/Au111Surface.dat

@@ -0,0 +1,192 @@
+#
+# File Name
+#
+
+System.CurrrentDirectory         ./    # default=./
+System.Name               Au111Surface
+level.of.stdout                   1    # default=1 (1-3)
+level.of.fileout                  0    # default=1 (1-3)
+
+#
+# Definition of Atomic Species
+#
+
+Species.Number       2
+<Definition.of.Atomic.Species
+   Au   Au7.0-s3p2d2   Au_CA19
+   H    H6.0-s2p1      H_CA19
+Definition.of.Atomic.Species>
+
+#
+# Atoms
+#
+
+Atoms.Number        4
+Atoms.SpeciesAndCoordinates.Unit   FRAC  # Ang|AU
+<Atoms.SpeciesAndCoordinates
+   1   Au    0.10000000000000    0.33333333333333    0.66666666666667     8.5     8.5     0.0     0.0     0.0     0.0     0     off
+   2   Au    0.11111111111111    0.66666666666667    0.33333333333333     8.5     8.5     0.0     0.0     0.0     0.0     0     off
+   3   Au    0.12222222222222    0.00000000000000    0.00000000000000     8.5     8.5     0.0     0.0     0.0     0.0     0     off
+   4   H    0.13333333333333    0.33333333333333    0.66666666666667     0.5     0.5     0.0     0.0     0.0     0.0     0     off
+Atoms.SpeciesAndCoordinates>
+   #5   H    0.14444444444444    0.66666666666667    0.33333333333333     0.5     0.5     0.0     0.0     0.0     0.0     0     off
+   #6   H    0.15555555555555    0.00000000000000    0.00000000000000     0.5     0.5     0.0     0.0     0.0     0.0     0     off
+
+Atoms.UnitVectors.Unit             Ang # Ang|AU
+<Atoms.UnitVectors
+  211.862723    0.000000    0.000000
+    0.000000    2.883086    0.000000
+    0.000000   -1.441543    2.496826
+Atoms.UnitVectors>
+
+#
+# SCF or Electronic System
+#
+
+scf.XcType                 LSDA-CA     # LDA|LSDA-CA|LSDA-PW
+scf.SpinPolarization        NC         # On|Off|NC
+scf.partialCoreCorrection   on         # On|Off
+scf.Hubbard.U		    off	       # On|Off , default=off 
+scf.Hubbard.Occupation     dual	       # onsite|full|dual , default=dual
+scf.SpinOrbit.Coupling      on         # On|Off , default=off
+scf.Constraint.NC.Spin      off        # On|Off , default=off
+scf.Constraint.NC.Spin.v    1.0        # default=0.0(ev)
+scf.ElectronicTemperature  100.0       # default=300 (K)
+scf.energycutoff           200.0       # default=150 (Ry)
+scf.maxIter                500         # default=40
+scf.EigenvalueSolver       Band        # Recursion|Cluster|Band
+scf.Kgrid                  1 12 12     # means n1 x n2 x n3
+scf.ProExpn.VNA             on         # on|off, default=on
+scf.Mixing.Type           Rmm-Diisk    # Simple|Rmm-Diis|Gr-Pulay
+scf.Init.Mixing.Weight     0.01        # default=0.30
+scf.Min.Mixing.Weight      0.001       # default=0.001
+scf.Max.Mixing.Weight      0.01        # default=0.40
+scf.Mixing.History          30         # default=5
+scf.Mixing.StartPulay       20         # default=6
+scf.criterion             1.0e-8       # default=1.0e-6 (Hartree)
+ESM.switch                 off         # off, on1=v|v|v, on2=m|v|m, on3=v|v|m, on4=on2+EF
+ESM.buffer.range           4.5         # default=10.0 (ang)
+ESM.wall.switch            on
+ESM.wall.position          6.7         # default=10.0 (ang)
+ESM.wall.height           100.0        # default=100.0 (eV)
+
+#
+# 1D FFT
+#
+
+1DFFT.NumGridK             900         # default=900
+1DFFT.NumGridR             900         # default=900
+1DFFT.EnergyCutoff        3600.0       # default=3DFFT.EnergyCutoff*3.0 (Ry)
+
+#
+# Orbital Optimization
+#
+
+orbitalOpt.Method           off        # Off|Unrestricted|Restricted
+orbitalOpt.InitCoes     Symmetrical    # Symmetrical|Free
+orbitalOpt.initPrefactor   0.1         # default=0.1
+orbitalOpt.scf.maxIter      15         # default=12
+orbitalOpt.MD.maxIter        5         # default=5
+orbitalOpt.per.MDIter       20         # default=1000000
+orbitalOpt.criterion      1.0e-4       # default=1.0e-4 (Hartree/borh)^2
+
+#
+# output of contracted orbitals
+#
+
+CntOrb.fileout               off       # on|off, default=off
+Num.CntOrb.Atoms              3        # default=1
+<Atoms.Cont.Orbitals
+ 1
+ 5
+ 9
+Atoms.Cont.Orbitals>
+
+#
+# SCF Order-N
+#
+
+orderN.HoppingRanges        9.2        # default=5.0 (Ang) 
+orderN.NumHoppings           3         # default=2
+orderN.KrylovH.order        600        # default=400
+orderN.recalc.EM            off
+
+#
+# restart using *.rst
+#
+
+scf.restart    off
+
+
+#
+# MD or Geometry Optimization
+#
+
+MD.Type                       Nomd     # Nomd|Constant_Energy_MD|Opt
+MD.maxIter                    1000     # default=1
+MD.Opt.DIIS.History           3        # default=3
+MD.Opt.StartDIIS              5        # default=5
+MD.TimeStep                   1        # default=0.5 (fs)
+MD.Opt.criterion            1.0e-4     # default=1.0e-4 (Hartree/bohr)
+
+#
+# Band dispersion
+#
+
+Band.dispersion              off       # on|off, default=off
+#Band.Nkpath                3
+#<Band.kpath
+#  135  0.0  0.500000  0.000000    0.0   0.000000  0.000000     M  G
+#  135  0.0  0.000000  0.000000    0.0   0.333333  0.333333     G  K
+#  135  0.0  0.333333  0.333333    0.0   0.500000  0.000000     K  M
+#Band.kpath>
+
+
+#
+# MO output
+#
+
+MO.fileout                       off   # on|off
+num.HOMOs                         3    # default=2
+num.LUMOs                         3    # default=2
+
+MO.Nkpoint                        1    # default=1
+<MO.kpoint
+  0.0  0.0  0.0
+MO.kpoint>
+
+#
+# DOS and PDOS
+#
+
+DosGauss.fileout	      off      # on|off , default=off
+DosGauss.Num.Mesh	     500       # default=200
+DosGauss.Width		     0.1       # default=0.2(eV)	
+
+Dos.fileout                   off      # on|off, default=off
+Dos.Erange              -10.0  10.0    # default = -20 20 
+Dos.Kgrid                12 12 12      # default = Kgrid1 Kgrid2 Kgrid3
+
+
+HS.fileout                    on       # on|off, default=off
+
+
+#
+# Spin textures; Mulliken population
+#
+
+Filename.scfout    Au111Surface.scfout # default: default
+Filename.outdata   Au111Surface_BD     # default: default
+Calc.Type          BandDispersion      # default: MulPOnly
+Energy.Range        -1.0  1.0          # eV; default: 0.0  0.0
+Band.Nkpath                2
+<Band.kpath
+  135  0.0  0.500000  0.000000    0.0   0.000000  0.000000     M  G
+  135  0.0  0.000000  0.000000    0.0  -0.500000  0.000000     G -M
+Band.kpath>
+Filename.atomMulP Au111Surface_BD.AMulPBand # default: default
+Filename.xyzdata   Au111Surface_BD_MC       # default: default
+Num.of.Extract.Atom        3                # default: 1
+Extract.Atom               1 2 3            # default: 1 2 ... (Num.of.Extract.Atom)
+MulP.Vec.Scale       0.1  0.1  0.1          # default: 1.0  1.0  1.0
+Data.Reduction             1                # default: 1

tests/openmx/Au111Surface.md

@@ -0,0 +1,12 @@
+4 
+   time=    0.000 (fs)  Energy= -317.65405 (Hartree) Cell_Vectors= 211.86272  0.00000  0.00000  0.00000  2.88309  0.00000  0.00000 -1.44154  2.49683  scf_conv=1
+  Au   21.18627 -0.00000  1.66455   0.00620  0.00000 -0.00012  0.00000  0.00000  0.00000   0.00297  0.00000 159.37740 172.68083
+  Au   23.54030  1.44154  0.83228   0.00359 -0.00000  0.00000  0.00000  0.00000  0.00000  -0.00169  0.00000 156.20488 170.35407
+  Au   25.89433  0.00000  0.00000   0.01489 -0.00000  0.00011  0.00000  0.00000  0.00000   0.03068  0.00000 170.57058 146.22112
+   H   28.24836 -0.00000  1.66455  -0.02466 -0.00000  0.00001  0.00000  0.00000  0.00000  -0.03197  0.00000 16.47069 -11.89518
+4 
+   time=    1.000 (fs)  Energy= -317.65405 (Hartree) Cell_Vectors= 211.86272  0.00000  0.00000  0.00000  2.88309  0.00000  0.00000 -1.44154  2.49683  scf_conv=0
+  Au   21.18627 -0.00000  1.66455  -0.12489 -0.00000 -0.00012  0.00000  0.00000  0.00000   0.00297  0.00000 159.37740 172.68083
+  Au   23.54030  1.44154  0.83228  -0.00359 -0.00000  0.00000  0.00000  0.00000  0.00000  -0.00169  0.00000 156.20488 170.35407
+  Au   25.89433  0.00000  0.00000   0.07945  0.00000  0.00011  0.00000  0.00000  0.00000   0.03068  0.00000 170.57058 146.22112
+   H   28.24836 -0.00000  1.66455   0.04904  0.00000  0.00001  0.00000  0.00000  0.00000  -0.03197  0.00000 16.47069 -11.89518