diff --git a/dpdata/openmx/omx.py b/dpdata/openmx/omx.py index aae9b5786..89f853687 100644 --- a/dpdata/openmx/omx.py +++ b/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) diff --git a/tests/openmx/Au111Surface.dat b/tests/openmx/Au111Surface.dat new file mode 100644 index 000000000..ec6c816b2 --- /dev/null +++ b/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 + + +# +# Atoms +# + +Atoms.Number 4 +Atoms.SpeciesAndCoordinates.Unit FRAC # Ang|AU + + #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 + + +# +# 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 + + +# +# 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 +# + + +# +# MO output +# + +MO.fileout off # on|off +num.HOMOs 3 # default=2 +num.LUMOs 3 # default=2 + +MO.Nkpoint 1 # default=1 + + +# +# 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 + +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 diff --git a/tests/openmx/Au111Surface.md b/tests/openmx/Au111Surface.md new file mode 100644 index 000000000..1ca5a00ff --- /dev/null +++ b/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 diff --git a/tests/openmx/Methane2.dat b/tests/openmx/Methane2.dat deleted file mode 100644 index f48bb5c31..000000000 --- a/tests/openmx/Methane2.dat +++ /dev/null @@ -1,68 +0,0 @@ -# -# File Name -# - -System.CurrrentDirectory ./ # default=./ -System.Name Methane2 -level.of.stdout 1 # default=1 (1-3) -level.of.fileout 1 # default=1 (0-2) - -# -# Definition of Atomic Species -# - -Species.Number 2 - - -# -# Atoms -# - -Atoms.Number 5 -Atoms.SpeciesAndCoordinates.Unit Ang # Ang|AU - -Atoms.UnitVectors.Unit Ang # Ang|AU -# - -# -# SCF or Electronic System -# - -scf.XcType GGA-PBE # LDA|LSDA-CA|LSDA-PW|GGA-PBE -scf.SpinPolarization off # On|Off|NC -scf.ElectronicTemperature 100.0 # default=300 (K) -scf.energycutoff 200.0 # default=150 (Ry) -scf.maxIter 1 # default=40 -scf.EigenvalueSolver cluster # DC|GDC|Cluster|Band -scf.Kgrid 1 1 1 # means n1 x n2 x n3 -scf.Mixing.Type rmm-diis # Simple|Rmm-Diis|Gr-Pulay|Kerker|Rmm-Diisk -scf.Init.Mixing.Weight 0.30 # default=0.30 -scf.Min.Mixing.Weight 0.001 # default=0.001 -scf.Max.Mixing.Weight 0.400 # default=0.40 -scf.Mixing.History 15 # default=5 -scf.Mixing.StartPulay 4 # default=6 -scf.criterion 1.0e-8 # default=1.0e-6 (Hartree) - -# -# MD or Geometry Optimization -# - -MD.Type opt # Nomd|Opt|DIIS|NVE|NVT_VS|NVT_NH -MD.Opt.DIIS.History 7 # default=7 -MD.Opt.StartDIIS 5 # default=5 -MD.maxIter 5 # default=1 -MD.TimeStep 1.0 # default=0.5 (fs) -MD.Opt.criterion 1.0e-4 # default=1.0e-4 (Hartree/bohr) diff --git a/tests/openmx/Methane2.md b/tests/openmx/Methane2.md deleted file mode 100644 index 25ddedb84..000000000 --- a/tests/openmx/Methane2.md +++ /dev/null @@ -1,35 +0,0 @@ -5 - time= 0.000 (fs) Energy= -8.15440 (Hartree) Cell_Vectors= 10.00000 0.00000 0.00000 0.00000 10.00000 0.00000 0.00000 0.00000 10.00000 - C 0.30000 0.00000 0.00000 -0.36382 0.22843 -0.00000 0.00000 0.00000 0.00000 0.17513 0.00000 0.00000 0.00000 0 - H -0.88998 -0.62931 0.00000 0.04918 0.01544 0.00000 0.00000 0.00000 0.00000 -0.07837 0.00000 0.00000 0.00000 0 - H 0.00000 0.62931 -0.88998 0.02120 -0.00206 -0.00338 0.00000 0.00000 0.00000 -0.01532 0.00000 0.00000 0.00000 0 - H 0.00000 0.62931 0.88998 0.02120 -0.00206 0.00338 0.00000 0.00000 0.00000 -0.01532 0.00000 0.00000 0.00000 0 - H 0.88998 -0.62931 0.00000 0.23151 -0.22432 -0.00000 0.00000 0.00000 0.00000 -0.06612 0.00000 0.00000 0.00000 0 -5 - time= 1.000 (fs) Energy= -8.22382 (Hartree) Cell_Vectors= 10.00000 0.00000 0.00000 0.00000 10.00000 0.00000 0.00000 0.00000 10.00000 - C 0.21037 0.05628 -0.00000 -0.12208 -0.03279 0.00000 0.00000 0.00000 0.00000 0.18366 0.00000 0.00000 0.00000 0 - H -0.87786 -0.62551 0.00000 0.03791 0.01857 0.00000 0.00000 0.00000 0.00000 -0.08412 0.00000 0.00000 0.00000 0 - H 0.00522 0.62881 -0.89081 0.00960 0.01866 -0.02578 0.00000 0.00000 0.00000 -0.03522 0.00000 0.00000 0.00000 0 - H 0.00522 0.62881 0.89081 0.00960 0.01866 0.02578 0.00000 0.00000 0.00000 -0.03522 0.00000 0.00000 0.00000 0 - H 0.94702 -0.68458 -0.00000 0.04642 -0.02939 -0.00000 0.00000 0.00000 0.00000 -0.02910 0.00000 0.00000 0.00000 0 -5 - time= 2.000 (fs) Energy= -8.24265 (Hartree) Cell_Vectors= 10.00000 0.00000 0.00000 0.00000 10.00000 0.00000 0.00000 0.00000 10.00000 - C 0.12898 0.03442 0.00000 -0.04237 -0.04096 0.00000 0.00000 0.00000 0.00000 0.18611 0.00000 0.00000 0.00000 0 - H -0.85259 -0.61313 0.00000 0.01505 0.00643 0.00000 0.00000 0.00000 0.00000 -0.06398 0.00000 0.00000 0.00000 0 - H 0.01162 0.64125 -0.90800 0.00781 0.01155 -0.01439 0.00000 0.00000 0.00000 -0.04343 0.00000 0.00000 0.00000 0 - H 0.01162 0.64125 0.90800 0.00781 0.01155 0.01439 0.00000 0.00000 0.00000 -0.04343 0.00000 0.00000 0.00000 0 - H 0.97796 -0.70417 -0.00000 0.00515 0.00555 -0.00000 0.00000 0.00000 0.00000 -0.03526 0.00000 0.00000 0.00000 0 -5 - time= 3.000 (fs) Energy= -8.24626 (Hartree) Cell_Vectors= 10.00000 0.00000 0.00000 0.00000 10.00000 0.00000 0.00000 0.00000 10.00000 - C 0.10073 0.00711 0.00000 -0.01636 -0.00597 0.00000 0.00000 0.00000 0.00000 0.18796 0.00000 0.00000 0.00000 0 - H -0.84256 -0.60884 0.00000 -0.00199 -0.00521 -0.00000 0.00000 0.00000 0.00000 -0.05693 0.00000 0.00000 0.00000 0 - H 0.01683 0.64895 -0.91760 0.00694 0.00262 -0.00333 0.00000 0.00000 0.00000 -0.04635 0.00000 0.00000 0.00000 0 - H 0.01683 0.64895 0.91760 0.00694 0.00262 0.00333 0.00000 0.00000 0.00000 -0.04635 0.00000 0.00000 0.00000 0 - H 0.98139 -0.70048 -0.00000 0.00184 0.00495 -0.00000 0.00000 0.00000 0.00000 -0.03833 0.00000 0.00000 0.00000 0 -5 - time= 4.000 (fs) Energy= -8.24682 (Hartree) Cell_Vectors= 10.00000 0.00000 0.00000 0.00000 10.00000 0.00000 0.00000 0.00000 10.00000 - C 0.08982 0.00313 0.00000 -0.00704 -0.00134 0.00000 0.00000 0.00000 0.00000 0.19038 0.00000 0.00000 0.00000 0 - H -0.84389 -0.61232 0.00000 -0.00554 -0.00705 0.00000 0.00000 0.00000 0.00000 -0.05558 0.00000 0.00000 0.00000 0 - H 0.02145 0.65069 -0.91982 0.00591 0.00121 -0.00163 0.00000 0.00000 0.00000 -0.04727 0.00000 0.00000 0.00000 0 - H 0.02145 0.65069 0.91982 0.00591 0.00121 0.00163 0.00000 0.00000 0.00000 -0.04727 0.00000 0.00000 0.00000 0 - H 0.98262 -0.69717 -0.00000 -0.00051 0.00564 -0.00000 0.00000 0.00000 0.00000 -0.04026 0.00000 0.00000 0.00000 0 diff --git a/tests/test_openmx_check_convergence.py b/tests/test_openmx_check_convergence.py deleted file mode 100644 index b19ad6e8d..000000000 --- a/tests/test_openmx_check_convergence.py +++ /dev/null @@ -1,64 +0,0 @@ -from __future__ import annotations - -import unittest - -import numpy as np -from context import dpdata - - -class TestOPENMXTRAJProps: - def test_atom_names(self): - self.assertEqual(self.system.data["atom_names"], ["C", "H"]) - - def test_atom_numbs(self): - self.assertEqual(self.system.data["atom_numbs"], [1, 4]) - - def test_atom_types(self): - for ii in range(0, 1): - self.assertEqual(self.system.data["atom_types"][ii], 0) - for ii in range(1, 5): - self.assertEqual(self.system.data["atom_types"][ii], 1) - - def test_cell(self): - ref = 10.0 * np.eye(3) - self.assertEqual(self.system.get_nframes(), 5) - for ff in range(self.system.get_nframes()): - for ii in range(3): - for jj in range(3): - self.assertEqual(self.system["cells"][ff][ii][jj], ref[ii][jj]) - - def test_coord(self): - with open("openmx/Methane2.md") as md_file: - lines = md_file.readlines() - lines = lines[-5:] - coords = [] - for line in lines: - parts = line.split() - for_line = [float(parts[1]), float(parts[2]), float(parts[3])] - coords.append(for_line) - coords = np.array(coords) - celll = 10.0 - ## Applying PBC ## - for ii in range(5): - for jj in range(3): - if coords[ii][jj] < 0: - coords[ii][jj] += celll - elif coords[ii][jj] >= celll: - coords[ii][jj] -= celll - self.assertAlmostEqual( - self.system["coords"][-1][ii][jj], coords[ii][jj] - ) - - -class TestOPENMXTraj(unittest.TestCase, TestOPENMXTRAJProps): - def setUp(self): - self.system = dpdata.System("openmx/Methane2", fmt="openmx/md") - - -class TestOPENMXLabeledTraj(unittest.TestCase, TestOPENMXTRAJProps): - def setUp(self): - self.system = dpdata.LabeledSystem("openmx/Methane2", fmt="openmx/md") - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/test_openmx_multiple_formats.py b/tests/test_openmx_multiple_formats.py new file mode 100644 index 000000000..e1892a57c --- /dev/null +++ b/tests/test_openmx_multiple_formats.py @@ -0,0 +1,66 @@ +from __future__ import annotations + +import unittest +import warnings + +import numpy as np +from context import dpdata + + +class TestOPENMXTRAJProps: + def test_atom_names(self): + self.assertEqual(self.system.data["atom_names"], ["Au", "H"]) + + def test_atom_numbs(self): + self.assertEqual(self.system.data["atom_numbs"], [3, 1]) + + def test_atom_types(self): + self.assertEqual(list(self.system.data["atom_types"]), [0, 0, 0, 1]) + + def test_cell(self): + cell = np.array( + [[211.86272, 0.0, 0.0], [0.0, 2.88309, 0.0], [0.0, -1.44154, 2.49683]] + ) + cells = np.array([cell, cell]) + self.assertEqual(self.system.get_nframes(), 2) + for ff in range(self.system.get_nframes()): + for ii in range(3): + for jj in range(3): + self.assertAlmostEqual( + self.system["cells"][ff][ii][jj], cells[ff][ii][jj] + ) + + def test_coord(self): + coord = np.array( + [ + [21.18627, 0.0, 1.66455], + [23.5403, 1.44154, 0.83228], + [25.89433, 0.0, 0.0], + [28.24836, 0.0, 1.66455], + ] + ) + coords = np.array([coord, coord]) + for ff in range(self.system.get_nframes()): + for ii in range(4): + for jj in range(3): + self.assertAlmostEqual( + self.system["coords"][ff][ii][jj], coords[ff][ii][jj] + ) + + +class TestOPENMXTraj(unittest.TestCase, TestOPENMXTRAJProps): + def setUp(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + self.system = dpdata.System("openmx/Au111Surface", fmt="openmx/md") + + +class TestOPENMXLabeledTraj(unittest.TestCase, TestOPENMXTRAJProps): + def setUp(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + self.system = dpdata.LabeledSystem("openmx/Au111Surface", fmt="openmx/md") + + +if __name__ == "__main__": + unittest.main()