diff --git a/docs/credits.rst b/docs/credits.rst new file mode 100644 index 000000000..a72b83e5a --- /dev/null +++ b/docs/credits.rst @@ -0,0 +1,4 @@ +Authors +======= + +.. git-shortlog-authors:: \ No newline at end of file diff --git a/docs/index.rst b/docs/index.rst index 64d9183af..e4b34f1db 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -13,6 +13,7 @@ Welcome to dpdata's documentation! cli formats api/api + credits .. mdinclude:: ../README.md diff --git a/dpdata/abacus/relax.py b/dpdata/abacus/relax.py index 3b334ab92..c9d7a669d 100644 --- a/dpdata/abacus/relax.py +++ b/dpdata/abacus/relax.py @@ -68,6 +68,7 @@ def get_coords_from_log(loglines,natoms): else: cells.append(cells[-1]) + coords[-1] = np.array(coords[-1]) if direct_coord: coords[-1] = coords[-1].dot(cells[-1]) diff --git a/dpdata/abacus/scf.py b/dpdata/abacus/scf.py index e6d76b635..45ca9a663 100644 --- a/dpdata/abacus/scf.py +++ b/dpdata/abacus/scf.py @@ -104,15 +104,11 @@ def get_energy(outlines): Etot = float(line.split()[1]) # in eV break if not Etot: - not_converge = False - for line in outlines: - if "convergence has NOT been achieved!" in line: - not_converge = True - raise RuntimeError("convergence has NOT been achieved in scf!") - break - if not not_converge: - raise RuntimeError("Final total energy cannot be found in output. Unknown problem.") - return Etot + raise RuntimeError("Final total energy cannot be found in output. Unknown problem.") + for line in outlines: + if "convergence has NOT been achieved!" in line: + return Etot,False + return Etot,True def get_force (outlines, natoms): force = [] @@ -156,7 +152,15 @@ def get_frame (fname): celldm, cell = get_cell(geometry_inlines) atom_names, natoms, types, coords = get_coords(celldm, cell, geometry_inlines, inlines) - energy = get_energy(outlines) + energy,converge = get_energy(outlines) + if not converge: + return {'atom_names':atom_names,\ + 'atom_numbs':natoms,\ + 'atom_types':types,\ + 'cells':[],\ + 'coords':[],\ + 'energies':[],\ + 'forces':[]} force = get_force (outlines, natoms) stress = get_stress(outlines) if stress is not None: diff --git a/dpdata/cp2k/output.py b/dpdata/cp2k/output.py index f91d86045..2f5e9cc61 100644 --- a/dpdata/cp2k/output.py +++ b/dpdata/cp2k/output.py @@ -104,8 +104,8 @@ def handle_single_log_frame(self, lines): # CONSERVED QUANTITY [hartree] = -0.279168013085E+04 energy_pattern_2 = re.compile(r' POTENTIAL ENERGY\[hartree\]\s+=\s+(?P\S+)') energy=None - cell_length_pattern = re.compile(r' INITIAL CELL LNTHS\[bohr\]\s+=\s+(?P\S+)\s+(?P\S+)\s+(?P\S+)') - cell_angle_pattern = re.compile(r' INITIAL CELL ANGLS\[deg\]\s+=\s+(?P\S+)\s+(?P\S+)\s+(?P\S+)') + cell_length_pattern = re.compile(r' (INITIAL ){0,1}CELL LNTHS\[bohr\]\s+=\s+(?P\S+)\s+(?P\S+)\s+(?P\S+)') + cell_angle_pattern = re.compile(r' (INITIAL ){0,1}CELL ANGLS\[deg\]\s+=\s+(?P\S+)\s+(?P\S+)\s+(?P\S+)') cell_A, cell_B, cell_C = (0,0,0,) cell_alpha, cell_beta, cell_gamma=(0,0,0,) cell_a_pattern = re.compile(r' CELL\| Vector a \[angstrom\]:\s+(?P\S+)\s+(?P\S+)\s+(?P\S+)') diff --git a/dpdata/driver.py b/dpdata/driver.py index 97583a101..670b03378 100644 --- a/dpdata/driver.py +++ b/dpdata/driver.py @@ -147,3 +147,74 @@ def label(self, data: dict) -> dict: labeled_data['energies'] += lb_data ['energies'] labeled_data['forces'] += lb_data ['forces'] return labeled_data + + +class Minimizer(ABC): + """The base class for a minimizer plugin. A minimizer can + minimize geometry. + """ + __MinimizerPlugin = Plugin() + + @staticmethod + def register(key: str) -> Callable: + """Register a minimizer plugin. Used as decorators. + + Parameter + --------- + key: str + key of the plugin. + + Returns + ------- + Callable + decorator of a class + + Examples + -------- + >>> @Minimizer.register("some_minimizer") + ... class SomeMinimizer(Minimizer): + ... pass + """ + return Minimizer.__MinimizerPlugin.register(key) + + @staticmethod + def get_minimizer(key: str) -> "Minimizer": + """Get a minimizer plugin. + + Parameter + --------- + key: str + key of the plugin. + + Returns + ------- + Minimizer + the specific minimizer class + + Raises + ------ + RuntimeError + if the requested minimizer is not implemented + """ + try: + return Minimizer.__MinimizerPlugin.plugins[key] + except KeyError as e: + raise RuntimeError('Unknown minimizer: ' + key) from e + + def __init__(self, *args, **kwargs) -> None: + """Setup the minimizer.""" + + @abstractmethod + def minimize(self, data: dict) -> dict: + """Minimize the geometry. + + Parameters + ---------- + data : dict + data with coordinates and atom types + + Returns + ------- + dict + labeled data with minimized coordinates, energies, and forces + """ diff --git a/dpdata/fhi_aims/output.py b/dpdata/fhi_aims/output.py index 4ee819c3b..1a1b2c579 100755 --- a/dpdata/fhi_aims/output.py +++ b/dpdata/fhi_aims/output.py @@ -1,5 +1,6 @@ import numpy as np import re +import warnings latt_patt="\|\s+([0-9]{1,}[.][0-9]*)\s+([0-9]{1,}[.][0-9]*)\s+([0-9]{1,}[.][0-9]*)" pos_patt_first="\|\s+[0-9]{1,}[:]\s\w+\s(\w+)(\s.*[-]?[0-9]{1,}[.][0-9]*)(\s+[-]?[0-9]{1,}[.][0-9]*)(\s+[-]?[0-9]{1,}[.][0-9]*)" @@ -63,7 +64,7 @@ def get_fhi_aims_block(fp) : return blk return blk -def get_frames (fname, md=True, begin = 0, step = 1) : +def get_frames (fname, md=True, begin = 0, step = 1, convergence_check=True) : fp = open(fname) blk = get_fhi_aims_block(fp) ret = get_info(blk, type_idx_zero = True) @@ -78,6 +79,7 @@ def get_frames (fname, md=True, begin = 0, step = 1) : all_virials = [] cc = 0 + rec_failed = [] while len(blk) > 0 : if debug: with open(str(cc),'w') as f: @@ -87,9 +89,9 @@ def get_frames (fname, md=True, begin = 0, step = 1) : coord, _cell, energy, force, virial, is_converge = analyze_block(blk, first_blk=True, md=md) else: coord, _cell, energy, force, virial, is_converge = analyze_block(blk, first_blk=False) - if is_converge : - if len(coord) == 0: - break + if len(coord) == 0: + break + if is_converge or not convergence_check: all_coords.append(coord) if _cell: @@ -101,9 +103,16 @@ def get_frames (fname, md=True, begin = 0, step = 1) : all_forces.append(force) if virial is not None : all_virials.append(virial) + if not is_converge: + rec_failed.append(cc+1) + blk = get_fhi_aims_block(fp) cc += 1 + 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 : diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py index 7136d3e63..e40dc400b 100644 --- a/dpdata/lammps/dump.py +++ b/dpdata/lammps/dump.py @@ -5,6 +5,11 @@ lib_path = os.path.dirname(os.path.realpath(__file__)) sys.path.append(lib_path) import lmp +import warnings +class UnwrapWarning(UserWarning): + pass +warnings.simplefilter('once', UnwrapWarning) + def _get_block (lines, key) : for idx in range(len(lines)) : @@ -59,16 +64,17 @@ def get_coordtype_and_scalefactor(keys): key_su = ['xsu','ysu','zsu'] #scaled and unfolded,sf = lattice parameter lmp_coor_type = [key_pc,key_uc,key_s,key_su] sf = [0,0,1,1] + uw = [0,1,0,1] # unwraped or not for k in range(4): if all(i in keys for i in lmp_coor_type[k]): - return lmp_coor_type[k],sf[k] + return lmp_coor_type[k], sf[k], uw[k] -def safe_get_posi(lines,cell,orig=np.zeros(3)) : +def safe_get_posi(lines,cell,orig=np.zeros(3), unwrap=False) : blk, head = _get_block(lines, 'ATOMS') keys = head.split() coord_tp_and_sf = get_coordtype_and_scalefactor(keys) assert coord_tp_and_sf is not None, 'Dump file does not contain atomic coordinates!' - coordtype, sf = coord_tp_and_sf + coordtype, sf, uw = coord_tp_and_sf id_idx = keys.index('id') - 2 xidx = keys.index(coordtype[0])-2 yidx = keys.index(coordtype[1])-2 @@ -81,8 +87,13 @@ def safe_get_posi(lines,cell,orig=np.zeros(3)) : posis.sort() posis = np.array(posis)[:,1:4] if not sf: - posis = (posis-orig)@np.linalg.inv(cell)# Convert to scaled coordinates for unscaled coordinates - return (posis%1)@cell # Convert scaled coordinates back to Cartesien coordinates + posis = (posis - orig) @ np.linalg.inv(cell) # Convert to scaled coordinates for unscaled coordinates + if uw and unwrap: + return posis @ cell # convert scaled coordinates back to Cartesien coordinates unwrap at the periodic boundaries + else: + if uw and not unwrap: + warnings.warn(message='Your dump file contains unwrapped coordinates, but you did not specify unwrapping (unwrap = True). The default is wrapping at periodic boundaries (unwrap = False).\n',category=UnwrapWarning) + return (posis % 1) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions def get_dumpbox(lines) : blk, h = _get_block(lines, 'BOX BOUNDS') @@ -145,9 +156,9 @@ def load_file(fname, begin = 0, step = 1) : cc += 1 if cc >= begin and (cc - begin) % step == 0 : buff.append(line) - -def system_data(lines, type_map = None, type_idx_zero = True) : + +def system_data(lines, type_map = None, type_idx_zero = True, unwrap=False) : array_lines = split_traj(lines) lines = array_lines[0] system = {} @@ -166,12 +177,12 @@ def system_data(lines, type_map = None, type_idx_zero = True) : system['cells'] = [np.array(cell)] natoms = sum(system['atom_numbs']) system['atom_types'] = get_atype(lines, type_idx_zero = type_idx_zero) - system['coords'] = [safe_get_posi(lines, cell, np.array(orig))] + system['coords'] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] for ii in range(1, len(array_lines)) : bounds, tilt = get_dumpbox(array_lines[ii]) orig, cell = dumpbox2box(bounds, tilt) system['cells'].append(cell) - system['coords'].append(safe_get_posi(array_lines[ii], cell, np.array(orig))) + system['coords'].append(safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap)) system['cells'] = np.array(system['cells']) system['coords'] = np.array(system['coords']) return system @@ -181,7 +192,7 @@ def split_traj(dump_lines) : marks = [] for idx,ii in enumerate(dump_lines) : if 'ITEM: TIMESTEP' in ii : - marks.append(idx) + marks.append(idx) if len(marks) == 0 : return None elif len(marks) == 1 : @@ -191,9 +202,9 @@ def split_traj(dump_lines) : ret = [] for ii in marks : ret.append(dump_lines[ii:ii+block_size]) - # for ii in range(len(marks)-1): + # for ii in range(len(marks)-1): # assert(marks[ii+1] - marks[ii] == block_size) - return ret + return ret return None diff --git a/dpdata/plugins/3dmol.py b/dpdata/plugins/3dmol.py new file mode 100644 index 000000000..c30893295 --- /dev/null +++ b/dpdata/plugins/3dmol.py @@ -0,0 +1,44 @@ +from typing import Tuple +import numpy as np + +from dpdata.format import Format +from dpdata.xyz.xyz import coord_to_xyz + + +@Format.register("3dmol") +class Py3DMolFormat(Format): + """3DMol format. + + To use this format, py3Dmol should be installed in advance. + """ + def to_system(self, + data: dict, + f_idx: int = 0, + size: Tuple[int] = (300,300), + style: dict = {"stick":{}, "sphere":{"radius":0.4}}, + **kwargs): + """Show 3D structure of a frame in jupyter. + + Parameters + ---------- + data : dict + system data + f_idx : int + frame index to show + size : tuple[int] + (width, height) of the widget + style : dict + style of 3DMol. Read 3DMol documentation for details. + + Examples + -------- + >>> system.to_3dmol() + """ + import py3Dmol + types = np.array(data['atom_names'])[data['atom_types']] + xyz = coord_to_xyz(data['coords'][f_idx], types) + viewer = py3Dmol.view(width=size[0], height=size[1]) + viewer.addModel(xyz, 'xyz') + viewer.setStyle(style.copy()) + viewer.zoomTo() + return viewer diff --git a/dpdata/plugins/__init__.py b/dpdata/plugins/__init__.py index 4887501ca..ca097fd75 100644 --- a/dpdata/plugins/__init__.py +++ b/dpdata/plugins/__init__.py @@ -14,6 +14,9 @@ importlib.import_module(module_name, PACKAGE_BASE) # https://setuptools.readthedocs.io/en/latest/userguide/entry_point.html -eps = metadata.entry_points().get('dpdata.plugins', []) +try: + eps = metadata.entry_points(group='dpdata.plugins') +except TypeError: + eps = metadata.entry_points().get('dpdata.plugins', []) for ep in eps: plugin = ep.load() diff --git a/dpdata/plugins/amber.py b/dpdata/plugins/amber.py index f9580f2fd..4fe41c1e4 100644 --- a/dpdata/plugins/amber.py +++ b/dpdata/plugins/amber.py @@ -5,7 +5,7 @@ import dpdata.amber.md import dpdata.amber.sqm from dpdata.format import Format -from dpdata.driver import Driver +from dpdata.driver import Driver, Minimizer @Format.register("amber/md") @@ -122,3 +122,21 @@ def label(self, data: dict) -> dict: ) from e labeled_system.append(dpdata.LabeledSystem(out_fn, fmt="sqm/out")) return labeled_system.data + + +@Minimizer.register("sqm") +class SQMMinimizer(Minimizer): + """SQM minimizer. + + Parameters + ---------- + maxcyc : int, default=1000 + maximun cycle to minimize + """ + def __init__(self, maxcyc=1000, *args, **kwargs) -> None: + assert maxcyc > 0, "maxcyc should be more than 0 to minimize" + self.driver = SQMDriver(maxcyc=maxcyc, **kwargs) + + def minimize(self, data: dict) -> dict: + # sqm has minimize feature + return self.driver.label(data) diff --git a/dpdata/plugins/ase.py b/dpdata/plugins/ase.py index f1745b472..f61d46f6c 100644 --- a/dpdata/plugins/ase.py +++ b/dpdata/plugins/ase.py @@ -1,10 +1,13 @@ -from dpdata.driver import Driver +from typing import TYPE_CHECKING, Type +from dpdata.driver import Driver, Minimizer from dpdata.format import Format import numpy as np import dpdata try: import ase.io from ase.calculators.calculator import PropertyNotImplementedError + if TYPE_CHECKING: + from ase.optimize.optimize import Optimizer except ImportError: pass @@ -46,6 +49,7 @@ def from_system(self, atoms: "ase.Atoms", **kwargs) -> dict: 'cells': np.array([cells]).astype('float32'), 'coords': np.array([coords]).astype('float32'), 'orig': np.zeros(3), + 'nopbc': not np.any(atoms.get_pbc()), } return info_dict @@ -203,3 +207,61 @@ def label(self, data: dict) -> dict: ls = dpdata.LabeledSystem(atoms, fmt="ase/structure", type_map=data['atom_names']) labeled_system.append(ls) return labeled_system.data + + +@Minimizer.register("ase") +class ASEMinimizer(Minimizer): + """ASE minimizer. + + Parameters + ---------- + driver : Driver + dpdata driver + optimizer : type, optional + ase optimizer class + fmax : float, optional, default=5e-3 + force convergence criterion + optimizer_kwargs : dict, optional + other parameters for optimizer + """ + def __init__(self, + driver: Driver, + optimizer: Type["Optimizer"] = None, + fmax: float = 5e-3, + optimizer_kwargs: dict = {}) -> None: + self.calculator = driver.ase_calculator + if optimizer is None: + from ase.optimize import LBFGS + self.optimizer = LBFGS + else: + self.optimizer = optimizer + self.optimizer_kwargs = { + "logfile": None, + **optimizer_kwargs.copy(), + } + self.fmax = fmax + + def minimize(self, data: dict) -> dict: + """Minimize the geometry. + + Parameters + ---------- + data : dict + data with coordinates and atom types + + Returns + ------- + dict + labeled data with minimized coordinates, energies, and forces + """ + system = dpdata.System(data=data) + # list[Atoms] + structures = system.to_ase_structure() + labeled_system = dpdata.LabeledSystem() + for atoms in structures: + atoms.calc = self.calculator + dyn = self.optimizer(atoms, **self.optimizer_kwargs) + dyn.run(fmax=self.fmax) + ls = dpdata.LabeledSystem(atoms, fmt="ase/structure", type_map=data['atom_names']) + labeled_system.append(ls) + return labeled_system.data \ No newline at end of file diff --git a/dpdata/plugins/fhi_aims.py b/dpdata/plugins/fhi_aims.py index 96e000c6b..b1805c4ef 100644 --- a/dpdata/plugins/fhi_aims.py +++ b/dpdata/plugins/fhi_aims.py @@ -4,7 +4,7 @@ @Format.register("fhi_aims/md") @Format.register("fhi_aims/output") class FhiMDFormat(Format): - def from_labeled_system(self, file_name, md=True, begin = 0, step = 1, **kwargs): + def from_labeled_system(self, file_name, md=True, begin = 0, step = 1, convergence_check=True, **kwargs): data = {} data['atom_names'], \ data['atom_numbs'], \ @@ -14,7 +14,7 @@ def from_labeled_system(self, file_name, md=True, begin = 0, step = 1, **kwargs) data['energies'], \ data['forces'], \ tmp_virial, \ - = dpdata.fhi_aims.output.get_frames(file_name, md = md, begin = begin, step = step) + = dpdata.fhi_aims.output.get_frames(file_name, md = md, begin = begin, step = step, convergence_check=convergence_check) if tmp_virial is not None : data['virials'] = tmp_virial return data diff --git a/dpdata/plugins/lammps.py b/dpdata/plugins/lammps.py index f1356bb1d..8f9962967 100644 --- a/dpdata/plugins/lammps.py +++ b/dpdata/plugins/lammps.py @@ -40,6 +40,7 @@ def from_system(self, type_map=None, begin=0, step=1, + unwrap=False, **kwargs): lines = dpdata.lammps.dump.load_file(file_name, begin=begin, step=step) - return dpdata.lammps.dump.system_data(lines, type_map) + return dpdata.lammps.dump.system_data(lines, type_map, unwrap=unwrap) diff --git a/dpdata/plugins/pwmat.py b/dpdata/plugins/pwmat.py index 7756e0c5c..3365806e5 100644 --- a/dpdata/plugins/pwmat.py +++ b/dpdata/plugins/pwmat.py @@ -11,7 +11,7 @@ @Format.register("pwmat/output") class PwmatOutputFormat(Format): @Format.post("rot_lower_triangular") - def from_labeled_system(self, file_name, begin=0, step=1, **kwargs): + def from_labeled_system(self, file_name, begin=0, step=1, convergence_check=True, **kwargs): data = {} data['atom_names'], \ data['atom_numbs'], \ @@ -21,7 +21,7 @@ def from_labeled_system(self, file_name, begin=0, step=1, **kwargs): data['energies'], \ data['forces'], \ tmp_virial \ - = dpdata.pwmat.movement.get_frames(file_name, begin=begin, step=step) + = dpdata.pwmat.movement.get_frames(file_name, begin=begin, step=step, convergence_check=convergence_check) if tmp_virial is not None: data['virials'] = tmp_virial # scale virial to the unit of eV diff --git a/dpdata/plugins/vasp.py b/dpdata/plugins/vasp.py index 0e0475151..07ec34f17 100644 --- a/dpdata/plugins/vasp.py +++ b/dpdata/plugins/vasp.py @@ -54,7 +54,7 @@ def to_system(self, data, frame_idx=0, **kwargs): @Format.register("vasp/outcar") class VASPOutcarFormat(Format): @Format.post("rot_lower_triangular") - def from_labeled_system(self, file_name, begin=0, step=1, **kwargs): + def from_labeled_system(self, file_name, begin=0, step=1, convergence_check=True, **kwargs): data = {} ml = kwargs.get("ml", False) data['atom_names'], \ @@ -65,7 +65,7 @@ def from_labeled_system(self, file_name, begin=0, step=1, **kwargs): data['energies'], \ data['forces'], \ tmp_virial, \ - = dpdata.vasp.outcar.get_frames(file_name, begin=begin, step=step, ml=ml) + = dpdata.vasp.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 diff --git a/dpdata/pwmat/movement.py b/dpdata/pwmat/movement.py index 121f38855..c39950f0a 100644 --- a/dpdata/pwmat/movement.py +++ b/dpdata/pwmat/movement.py @@ -1,5 +1,6 @@ import numpy as np from ..periodic_table import ELEMENTS +import warnings def system_info (lines, type_idx_zero = False) : atom_names = [] @@ -49,7 +50,7 @@ def get_movement_block(fp) : return blk # we assume that the force is printed ... -def get_frames (fname, begin = 0, step = 1) : +def get_frames (fname, begin = 0, step = 1, convergence_check=True) : fp = open(fname) blk = get_movement_block(fp) @@ -64,20 +65,28 @@ def get_frames (fname, begin = 0, step = 1) : 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) - if is_converge : - if len(coord) == 0: - break + 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_movement_block(fp) cc += 1 + + 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 diff --git a/dpdata/system.py b/dpdata/system.py index 3c9ccec8c..a2226d030 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -6,7 +6,7 @@ import dpdata.md.pbc from copy import deepcopy from enum import Enum, unique -from typing import Any, Tuple +from typing import Any, Tuple, Union from monty.json import MSONable from monty.serialization import loadfn,dumpfn from dpdata.periodic_table import Element @@ -17,7 +17,7 @@ import dpdata.plugins from dpdata.plugin import Plugin from dpdata.format import Format -from dpdata.driver import Driver +from dpdata.driver import Driver, Minimizer from dpdata.utils import ( elements_index_map, @@ -123,7 +123,7 @@ def check(self, system: "System"): elif isinstance(data, list): if len(shape) and shape[0] != len(data): raise DataError("Length of %s is %d, but expected %d" % (self.name, - len(shape), shape[0])) + len(data), shape[0])) else: raise RuntimeError("Unsupported type to check shape") elif self.required: @@ -176,6 +176,7 @@ def __init__ (self, begin = 0, step = 1, data = None, + convergence_check = True, **kwargs) : """ Constructor @@ -192,14 +193,49 @@ def __init__ (self, - ``deepmd/raw``: deepmd-kit raw - ``deepmd/npy``: deepmd-kit compressed format (numpy binary) - ``vasp/poscar``: vasp POSCAR + - ``vasp/contcar``: vasp contcar + - ``vasp/string``: vasp string + - ``vasp/outcar``: vasp outcar + - ``vasp/xml``: vasp xml - ``qe/cp/traj``: Quantum Espresso CP trajectory files. should have: file_name+'.in' and file_name+'.pos' - ``qe/pw/scf``: Quantum Espresso PW single point calculations. Both input and output files are required. If file_name is a string, it denotes the output file name. Input file name is obtained by replacing 'out' by 'in' from file_name. Or file_name is a list, with the first element being the input file name and the second element being the output filename. - ``abacus/scf``: ABACUS pw/lcao scf. The directory containing INPUT file is required. - ``abacus/md``: ABACUS pw/lcao MD. The directory containing INPUT file is required. - - ``abacus/relax``: ABACUS pw/lcao relax or cell-relax. The directory containing INPUT file is required. + - ``abacus/relax``: ABACUS pw/lcao relax or cell-relax. The directory containing INPUT file is required. + - ``abacus/stru``: abacus stru + - ``abacus/lcao/scf``: abacus lcao scf + - ``abacus/pw/scf``: abacus pw scf + - ``abacus/lcao/md``: abacus lcao md + - ``abacus/pw/md``: abacus pw md + - ``abacus/lcao/relax``: abacus lcao relax + - ``abacus/pw/relax``: abacus pw relax - ``siesta/output``: siesta SCF output file - ``siesta/aimd_output``: siesta aimd output file - ``pwmat/atom.config``: pwmat atom.config + - ``pwmat/movement``: pwmat movement + - ``pwmat/output``: pwmat output + - ``pwmat/mlmd``: pwmat mlmd + - ``pwmat/final.config``: pwmat final.config + - ``quip/gap/xyz_file``: quip gap xyz_file + - ``quip/gap/xyz``: quip gap xyz + - ``fhi_aims/output``: fhi_aims output + - ``fhi_aims/md``: fhi_aims md + - ``fhi_aims/scf``: fhi_aims scf + - ``pymatgen/structure``: pymatgen structure + - ``pymatgen/molecule``: pymatgen molecule + - ``pymatgen/computedstructureentry``: pymatgen computedstructureentry + - ``amber/md``: amber md + - ``sqm/out``: sqm out + - ``sqm/in``: sqm in + - ``ase/structure``: ase structure + - ``gaussian/log``: gaussian log + - ``gaussian/md``: gaussian md + - ``gaussian/gjf``: gaussian gjf + - ``deepmd/comp``: deepmd comp + - ``deepmd/hdf5``: deepmd hdf5 + - ``gromacs/gro``: gromacs gro + - ``cp2k/aimd_output``: cp2k aimd_output + - ``cp2k/output``: cp2k output type_map : list of str Needed by formats lammps/lmp and lammps/dump. Maps atom type to name. The atom with type `ii` is mapped to `type_map[ii]`. If not provided the atom names are assigned to `'Type_1'`, `'Type_2'`, `'Type_3'`... @@ -208,7 +244,9 @@ def __init__ (self, step : int The number of skipped frames when loading MD trajectory. data : dict - The raw data of System class. + The raw data of System class. + convergence_check : boolean + Whether to request a convergence check. """ self.data = {} self.data['atom_numbs'] = [] @@ -224,7 +262,7 @@ def __init__ (self, return if file_name is None : return - self.from_fmt(file_name, fmt, type_map=type_map, begin= begin, step=step, **kwargs) + self.from_fmt(file_name, fmt, type_map=type_map, begin= begin, step=step, convergence_check=convergence_check, **kwargs) if type_map is not None: self.apply_type_map(type_map) @@ -331,7 +369,7 @@ def dump(self,filename,indent=4): dumpfn(self.as_dict(),filename,indent=indent) - def map_atom_types(self,type_map=None): + def map_atom_types(self, type_map=None) -> np.ndarray: """ Map the atom types of the system @@ -346,7 +384,7 @@ def map_atom_types(self,type_map=None): Returns ------- - new_atom_types : list + new_atom_types : np.ndarray The mapped atom types """ if isinstance(type_map,dict) or type_map is None: @@ -366,7 +404,7 @@ def map_atom_types(self,type_map=None): atom_types_list=[] for name, numb in zip(self.get_atom_names(), self.get_atom_numbs()): atom_types_list.extend([name]*numb) - new_atom_types=np.array([type_map[ii] for ii in atom_types_list],dtype=np.int) + new_atom_types = np.array([type_map[ii] for ii in atom_types_list], dtype=int) return new_atom_types @@ -831,6 +869,28 @@ def predict(self, *args: Any, driver: str="dp", **kwargs: Any) -> "LabeledSystem data = driver.label(self.data.copy()) return LabeledSystem(data=data) + def minimize(self, *args: Any, minimizer: Union[str, Minimizer], **kwargs: Any) -> "LabeledSystem": + """Minimize the geometry. + + Parameters + ---------- + *args : iterable + Arguments passing to the minimizer + minimizer : str or Minimizer + The assigned minimizer + **kwargs : dict + Other arguments passing to the minimizer + + Returns + ------- + labeled_sys : LabeledSystem + A new labeled system. + """ + if not isinstance(minimizer, Minimizer): + minimizer = Minimizer.get_minimizer(minimizer)(*args, **kwargs) + data = minimizer.minimize(self.data.copy()) + return LabeledSystem(data=data) + def pick_atom_idx(self, idx, nopbc=None): """Pick atom index @@ -1270,10 +1330,43 @@ def predict(self, *args: Any, driver="dp", **kwargs: Any) -> "MultiSystems": """ if not isinstance(driver, Driver): driver = Driver.get_driver(driver)(*args, **kwargs) - new_multisystems = dpdata.MultiSystems() + new_multisystems = dpdata.MultiSystems(type_map=self.atom_names) for ss in self: new_multisystems.append(ss.predict(*args, driver=driver, **kwargs)) return new_multisystems + + def minimize(self, *args: Any, minimizer: Union[str, Minimizer], **kwargs: Any) -> "MultiSystems": + """ + Minimize geometry by a minimizer. + + Parameters + ---------- + *args : iterable + Arguments passing to the minimizer + minimizer : str or Minimizer + The assigned minimizer + **kwargs : dict + Other arguments passing to the minimizer + + Returns + ------- + MultiSystems + A new labeled MultiSystems. + + Examples + -------- + Minimize a system using ASE BFGS along with a DP driver: + >>> from dpdata.driver import Driver + >>> from ase.optimize import BFGS + >>> driver = driver.get_driver("dp")("some_model.pb") + >>> some_system.minimize(minimizer="ase", driver=driver, optimizer=BFGS, fmax=1e-5) + """ + if not isinstance(minimizer, Minimizer): + minimizer = Minimizer.get_minimizer(minimizer)(*args, **kwargs) + new_multisystems = dpdata.MultiSystems(type_map=self.atom_names) + for ss in self: + new_multisystems.append(ss.minimize(*args, minimizer=minimizer, **kwargs)) + return new_multisystems def pick_atom_idx(self, idx, nopbc=None): """Pick atom index diff --git a/dpdata/vasp/outcar.py b/dpdata/vasp/outcar.py index 1efdea2d8..3e32a1461 100644 --- a/dpdata/vasp/outcar.py +++ b/dpdata/vasp/outcar.py @@ -1,5 +1,6 @@ import numpy as np import re +import warnings def system_info(lines, type_idx_zero = False): atom_names = [] @@ -52,7 +53,7 @@ def get_outcar_block(fp, ml = False): return blk # we assume that the force is printed ... -def get_frames(fname, begin = 0, step = 1, ml = False): +def get_frames(fname, begin = 0, step = 1, ml = False, convergence_check=True): fp = open(fname) blk = get_outcar_block(fp) @@ -66,22 +67,29 @@ def get_frames(fname, begin = 0, step = 1, ml = False): 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 is_converge : - if len(coord) == 0: - break + 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 - + + 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 : @@ -101,8 +109,8 @@ def analyze_block(lines, ntot, nelm, ml = False): #select different searching tokens based on the ml label energy_token = ['free energy TOTEN', 'free energy ML TOTEN'] energy_index = [4, 5] - viral_token = ['FORCE on cell =-STRESS in cart. coord. units', 'ML FORCE'] - viral_index = [14, 4] + 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) @@ -121,8 +129,12 @@ def analyze_block(lines, ntot, nelm, ml = False): tmp_l = lines[idx+cell_index[ml_index]+dd] cell.append([float(ss) for ss in tmp_l.replace('-',' -').split()[0:3]]) - elif viral_token[ml_index] in ii: - tmp_v = [float(ss) for ss in lines[idx+viral_index[ml_index]].split()[2:8]] + 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] diff --git a/setup.py b/setup.py index d04124ab5..4407c625d 100644 --- a/setup.py +++ b/setup.py @@ -54,7 +54,15 @@ 'ase': ['ase'], 'amber': ['parmed'], 'pymatgen': ['pymatgen'], - 'docs': ['sphinx', 'recommonmark', 'sphinx_rtd_theme>=1.0.0rc1', 'numpydoc', 'm2r2', 'deepmodeling-sphinx', 'sphinx-argparse'], + 'docs': [ + 'sphinx', + 'recommonmark', + 'sphinx_rtd_theme>=1.0.0rc1', + 'numpydoc', + 'm2r2', + 'deepmodeling-sphinx>=0.1.1', + 'sphinx-argparse', + ], }, entry_points={"console_scripts": ["dpdata = dpdata.cli:dpdata_cli"]}, ) diff --git a/tests/abacus.scf/INPUT.fail b/tests/abacus.scf/INPUT.fail new file mode 100644 index 000000000..d9fcc9378 --- /dev/null +++ b/tests/abacus.scf/INPUT.fail @@ -0,0 +1,18 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix ch4fail +stru_file STRU.ch4 #the filename of file containing atom positions +kpoint_file KPT.ch4 #the name of file containing k points +pseudo_dir ./ +ntype 2 +nbands 8 +#Parameters (Accuracy) +ecutwfc 100 +symmetry 1 +scf_nmax 50 +smearing_method gauss #type of smearing: gauss; fd; fixed; mp; mp2; mv +smearing_sigma 0.01 +mixing_type plain +mixing_beta 0.5 +cal_force 1 +cal_stress 1 diff --git a/tests/abacus.scf/INPUT b/tests/abacus.scf/INPUT.ok similarity index 100% rename from tests/abacus.scf/INPUT rename to tests/abacus.scf/INPUT.ok diff --git a/tests/abacus.scf/OUT.ch4fail/running_scf.log b/tests/abacus.scf/OUT.ch4fail/running_scf.log new file mode 100644 index 000000000..745ca3b2a --- /dev/null +++ b/tests/abacus.scf/OUT.ch4fail/running_scf.log @@ -0,0 +1,598 @@ + + WELCOME TO ABACUS + + 'Atomic-orbital Based Ab-initio Computation at UStc' + + Website: http://abacus.ustc.edu.cn/ + + Version: Parallel, v2.1.0 + Processor Number is 1 + Start Time is Fri May 7 16:18:50 2021 + + ------------------------------------------------------------------------------------ + + READING GENERAL INFORMATION + global_out_dir = OUT.ch4/ + global_in_card = INPUT + pseudo_dir = ./ + pseudo_type = auto + DRANK = 1 + DSIZE = 1 + DCOLOR = 1 + GRANK = 1 + GSIZE = 1 + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Reading atom information in unitcell: | + | From the input file and the structure file we know the number of | + | different elments in this unitcell, then we list the detail | + | information for each element, especially the zeta and polar atomic | + | orbital number for each element. The total atom number is counted. | + | We calculate the nearest atom distance for each atom and show the | + | Cartesian and Direct coordinates for each atom. We list the file | + | address for atomic orbitals. The volume and the lattice vectors | + | in real and reciprocal space is also shown. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + READING UNITCELL INFORMATION + ntype = 2 + atom label for species 1 = C + atom label for species 2 = H + lattice constant (Bohr) = 10 + lattice constant (Angstrom) = 5.29177 + + READING ATOM TYPE 1 + atom label = C + start magnetization = FALSE + L=0, number of zeta = 1 + L=1, number of zeta = 1 + L=2, number of zeta = 1 + number of atom for this type = 1 + + READING ATOM TYPE 2 + atom label = H + start magnetization = FALSE + L=0, number of zeta = 1 + L=1, number of zeta = 1 + L=2, number of zeta = 1 + number of atom for this type = 4 + + TOTAL ATOM NUMBER = 5 + + CARTESIAN COORDINATES ( UNIT = 10 Bohr ). + atom x y z mag + tauc_C1 0.981274803 0.861285385001 0.838442496 0 + tauc_H1 0.0235572019992 0.758025625 0.663513359999 0 + tauc_H2 0.78075702 0.889445934999 0.837363467999 0 + tauc_H3 0.0640916129996 0.0434389050006 0.840995502 0 + tauc_H4 0.0393212140007 0.756530859 0.00960920699981 0 + + + Volume (Bohr^3) = 1000 + Volume (A^3) = 148.184534296 + + Lattice vectors: (Cartesian coordinate: in unit of a_0) + +1 +0 +0 + +0 +1 +0 + +0 +0 +1 + Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0) + +1 +0 +0 + +0 +1 +0 + +0 -0 +1 + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Reading pseudopotentials files: | + | The pseudopotential file is in UPF format. The 'NC' indicates that | + | the type of pseudopotential is 'norm conserving'. Functional of | + | exchange and correlation is decided by 4 given parameters in UPF | + | file. We also read in the 'core correction' if there exists. | + | Also we can read the valence electrons number and the maximal | + | angular momentum used in this pseudopotential. We also read in the | + | trail wave function, trail atomic density and local-pseudopotential| + | on logrithmic grid. The non-local pseudopotential projector is also| + | read in if there is any. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + PAO radial cut off (Bohr) = 15 + + Read in pseudopotential file is C_ONCV_PBE-1.0.upf + pseudopotential type = NC + functional Ex = PBE + functional Ec = + functional GCEx = + functional GCEc = + nonlocal core correction = 0 + valence electrons = 4 + lmax = 1 + number of zeta = 0 + number of projectors = 4 + L of projector = 0 + L of projector = 0 + L of projector = 1 + L of projector = 1 + PAO radial cut off (Bohr) = 15 + + Read in pseudopotential file is H_ONCV_PBE-1.0.upf + pseudopotential type = NC + functional Ex = PBE + functional Ec = + functional GCEx = + functional GCEc = + nonlocal core correction = 0 + valence electrons = 1 + lmax = 0 + number of zeta = 0 + number of projectors = 2 + L of projector = 0 + L of projector = 0 + initial pseudo atomic orbital number = 0 + NLOCAL = 45 + + SETUP THE ELECTRONS NUMBER + electron number of element C = 4 + total electron number of element C = 4 + electron number of element H = 1 + total electron number of element H = 4 + occupied bands = 4 + NBANDS = 8 + DONE : SETUP UNITCELL Time : 0.0506949424744 (SEC) + + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Doing symmetry analysis: | + | We calculate the norm of 3 vectors and the angles between them, | + | the type of Bravais lattice is given. We can judge if the unticell | + | is a primitive cell. Finally we give the point group operation for | + | this unitcell. We we use the point group operations to do symmetry | + | analysis on given k-point mesh and the charge density. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + LATTICE VECTORS: (CARTESIAN COORDINATE: IN UNIT OF A0) + +1 +0 +0 + +0 +1 +0 + +0 +0 +1 + right hand lattice = 1 + NORM_A = 1 + NORM_B = 1 + NORM_C = 1 + ALPHA (DEGREE) = 90 + BETA (DEGREE) = 90 + GAMMA (DEGREE) = 90 + BRAVAIS TYPE = 1 + BRAVAIS LATTICE NAME = 01. Cubic P (simple) + STANDARD LATTICE VECTORS: (CARTESIAN COORDINATE: IN UNIT OF A0) + +1 +0 +0 + +0 +1 +0 + +0 +0 +1 + IBRAV = 1 + BRAVAIS = SIMPLE CUBIC + LATTICE CONSTANT A = 4.35889894354 + ibrav = 1 + ROTATION MATRICES = 48 + PURE POINT GROUP OPERATIONS = 1 + SPACE GROUP OPERATIONS = 1 + POINT GROUP = C_1 + DONE : SYMMETRY Time : 0.103345155716 (SEC) + + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Setup K-points | + | We setup the k-points according to input parameters. | + | The reduced k-points are set according to symmetry operations. | + | We treat the spin as another set of k-points. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + + SETUP K-POINTS + nspin = 1 + Input type of k points = Monkhorst-Pack(Gamma) + nkstot = 1 + nkstot_ibz = 1 + IBZ DirectX DirectY DirectZ Weight ibz2bz + 1 0 0 0 1 0 + nkstot now = 1 + + KPOINTS DIRECT_X DIRECT_Y DIRECT_Z WEIGHT + 1 0 0 0 1 + + k-point number in this process = 1 + minimum distributed K point number = 1 + + KPOINTS CARTESIAN_X CARTESIAN_Y CARTESIAN_Z WEIGHT + 1 0 0 0 2 + + KPOINTS DIRECT_X DIRECT_Y DIRECT_Z WEIGHT + 1 0 0 0 2 + DONE : INIT K-POINTS Time : 0.103604078293 (SEC) + + + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + | | + | Setup plane waves: | + | Use the energy cutoff and the lattice vectors to generate the | + | dimensions of FFT grid. The number of FFT grid on each processor | + | is 'nrxx'. The number of plane wave basis in reciprocal space is | + | different for charege/potential and wave functions. We also set | + | the 'sticks' for the parallel of FFT. | + | | + <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + + SETUP THE PLANE WAVE BASIS + energy cutoff for wavefunc (unit:Ry) = 100 + [fft grid for wave functions] = 64, 64, 64 + [fft grid for charge/potential] = 64, 64, 64 + [fft grid division] = 1, 1, 1 + [big fft grid for charge/potential] = 64, 64, 64 + nbxx = 262144 + nrxx = 262144 + + SETUP PLANE WAVES FOR CHARGE/POTENTIAL + number of plane waves = 135043 + number of sticks = 3181 + + SETUP PLANE WAVES FOR WAVE FUNCTIONS + number of plane waves = 16879 + number of sticks = 793 + + PARALLEL PW FOR CHARGE/POTENTIAL + PROC COLUMNS(POT) PW + 1 3181 135043 + --------------- sum ------------------- + 1 3181 135043 + + PARALLEL PW FOR WAVE FUNCTIONS + PROC COLUMNS(W) PW + 1 793 16879 + --------------- sum ------------------- + 1 793 16879 + + SETUP COORDINATES OF PLANE WAVES + number of total plane waves = 135043 + + SETUP COORDINATES OF PLANE WAVES + number of |g| = 847 + max |g| = 1013 + min |g| = 0 + DONE : INIT PLANEWAVE Time : 0.294428110123 (SEC) + + npwx = 16879 + + SETUP NONLOCAL PSEUDOPOTENTIALS IN PLANE WAVE BASIS + C non-local projectors: + projector 1 L=0 + projector 2 L=0 + projector 3 L=1 + projector 4 L=1 + H non-local projectors: + projector 1 L=0 + projector 2 L=0 + TOTAL NUMBER OF NONLOCAL PROJECTORS = 16 + DONE : LOCAL POTENTIAL Time : 0.339339256287 (SEC) + + + Init Non-Local PseudoPotential table : + Init Non-Local-Pseudopotential done. + DONE : NON-LOCAL POTENTIAL Time : 0.352857112885 (SEC) + + start_pot = atomic + DONE : INIT POTENTIAL Time : 0.642299 (SEC) + + + Make real space PAO into reciprocal space. + max mesh points in Pseudopotential = 601 + dq(describe PAO in reciprocal space) = 0.01 + max q = 1206 + + number of pseudo atomic orbitals for C is 0 + + number of pseudo atomic orbitals for H is 0 + DONE : INIT BASIS Time : 0.859031 (SEC) + + ------------------------------------------- + ------------------------------------------- + + PW ALGORITHM --------------- ION= 1 ELEC= 1-------------------------------- + K-point CG iter num Time(Sec) + 1 8.000000 1.417983 + + Density error is 0.705428908048 + Error Threshold = 0.010000000000 + + Energy Rydberg eV + E_KohnSham -16.017544057872 -217.929867153100 + E_Harris -16.408058532460 -223.243089158968 + E_Fermi -0.626082258799 -8.518286136372 + + PW ALGORITHM --------------- ION= 1 ELEC= 2-------------------------------- + K-point CG iter num Time(Sec) + 1 3.000000 0.623305 + + Density error is 0.037193130746 + Error Threshold = 0.008817861351 + + Energy Rydberg eV + E_KohnSham -16.144996950340 -219.663952717248 + E_Harris -16.158845333973 -219.852369642740 + E_Fermi -0.425221736464 -5.785438529359 + + PW ALGORITHM --------------- ION= 1 ELEC= 3-------------------------------- + K-point CG iter num Time(Sec) + 1 3.750000 0.765803 + + Density error is 0.008021287795 + Error Threshold = 0.000464914134 + + Energy Rydberg eV + E_KohnSham -16.143766641525 -219.647213507066 + E_Harris -16.146720471833 -219.687402430176 + E_Fermi -0.423547802407 -5.762663488108 + + PW ALGORITHM --------------- ION= 1 ELEC= 4-------------------------------- + K-point CG iter num Time(Sec) + 1 3.250000 0.691389 + + Density error is 0.001345350159 + Error Threshold = 0.000100266097 + + Energy Rydberg eV + E_KohnSham -16.143852570867 -219.648382635743 + E_Harris -16.144348491782 -219.655129985945 + E_Fermi -0.417940644498 -5.686374190966 + + PW ALGORITHM --------------- ION= 1 ELEC= 5-------------------------------- + K-point CG iter num Time(Sec) + 1 3.500000 0.724421 + + Density error is 0.000193438710 + Error Threshold = 0.000016816877 + + Energy Rydberg eV + E_KohnSham -16.143964486009 -219.649905319365 + E_Harris -16.144043239940 -219.650976821572 + E_Fermi -0.412189586104 -5.608127027270 + + PW ALGORITHM --------------- ION= 1 ELEC= 6-------------------------------- + K-point CG iter num Time(Sec) + 1 4.750000 0.929892 + + Density error is 0.000033144966 + Error Threshold = 0.000002417984 + + Energy Rydberg eV + E_KohnSham -16.143962062113 -219.649872340569 + E_Harris -16.143973121330 -219.650022808940 + E_Fermi -0.408955040431 -5.564118775681 + + PW ALGORITHM --------------- ION= 1 ELEC= 7-------------------------------- + K-point CG iter num Time(Sec) + 1 2.625000 0.571903 + + Density error is 0.000005013879 + Error Threshold = 0.000000414312 + + Energy Rydberg eV + E_KohnSham -16.143964909086 -219.649911075624 + E_Harris -16.143965842245 -219.649923771901 + E_Fermi -0.407628831587 -5.546074778663 + + PW ALGORITHM --------------- ION= 1 ELEC= 8-------------------------------- + K-point CG iter num Time(Sec) + 1 3.625000 0.729712 + + Density error is 0.000001011179 + Error Threshold = 0.000000062673 + + Energy Rydberg eV + E_KohnSham -16.143964999058 -219.649912299760 + E_Harris -16.143965129574 -219.649914075510 + E_Fermi -0.407392472263 -5.542858945081 + + PW ALGORITHM --------------- ION= 1 ELEC= 9-------------------------------- + K-point CG iter num Time(Sec) + 1 3.500000 0.698617 + + Density error is 0.000000165426 + Error Threshold = 0.000000012640 + + Energy Rydberg eV + E_KohnSham -16.143965114417 -219.649913869298 + E_Harris -16.143965110616 -219.649913817582 + E_Fermi -0.407267643914 -5.541160568271 + + PW ALGORITHM --------------- ION= 1 ELEC= 10-------------------------------- + K-point CG iter num Time(Sec) + 1 3.000000 0.629669 + + Density error is 0.000000034854 + Error Threshold = 0.000000002068 + + Energy Rydberg eV + E_KohnSham -16.143965129640 -219.649914076412 + E_Harris -16.143965128874 -219.649914065985 + E_Fermi -0.407195695883 -5.540181665084 + + PW ALGORITHM --------------- ION= 1 ELEC= 11-------------------------------- + K-point CG iter num Time(Sec) + 1 4.000000 0.813535 + + Density error is 0.000000008107 + Error Threshold = 0.000000000436 + + Energy Rydberg eV + E_KohnSham -16.143965128688 -219.649914063463 + E_Harris -16.143965130414 -219.649914086943 + E_Fermi -0.407157690601 -5.539664576698 + + PW ALGORITHM --------------- ION= 1 ELEC= 12-------------------------------- + K-point CG iter num Time(Sec) + 1 3.125000 0.674851 + + Density error is 0.000000001913 + Error Threshold = 0.000000000101 + + Energy Rydberg eV + E_KohnSham -16.143965127817 -219.649914051607 + E_Harris -16.143965128941 -219.649914066906 + E_Fermi -0.407141275858 -5.539441242659 + + PW ALGORITHM --------------- ION= 1 ELEC= 13-------------------------------- + K-point CG iter num Time(Sec) + 1 4.000000 0.804644 + + Density error is 0.000000000482 + Error Threshold = 0.000000000024 + + Energy Rydberg eV + E_KohnSham -16.143965127167 -219.649914042766 + E_Harris -16.143965127871 -219.649914052349 + E_band -5.963675568058 -81.139968748982 + E_one_elec -25.076858425354 -341.188162524121 + E_Hartree +13.701412040584 +186.417274397756 + E_xc -6.404563484424 -87.138556590897 + E_Ewald +1.636044742026 +22.259530674496 + E_demet -0.000000000000 -0.000000000000 + E_descf +0.000000000000 +0.000000000000 + E_efield +0.000000000000 +0.000000000000 + E_Fermi -0.407135332076 -5.539360373357 + convergence has NOT been achieved! + + STATE ENERGY(eV) AND OCCUPATIONS. 1/1 kpoint (Cartesian) = 0.00000 0.00000 0.00000 (16879 pws) + [spin1_state] 1 -15.903159 2.000000 + [spin1_state] 2 -8.412500 2.000000 + [spin1_state] 3 -8.255893 2.000000 + [spin1_state] 4 -7.998432 2.000000 + [spin1_state] 5 -0.514947 0.000000 + [spin1_state] 6 2.727369 0.000000 + [spin1_state] 7 3.067094 0.000000 + [spin1_state] 8 4.824439 0.000000 + + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (Ry/Bohr) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + C1 -0.006028 -0.043357 +0.003245 + H1 +0.011378 +0.004647 +0.013137 + H2 -0.026755 -0.000770 +0.009545 + H3 +0.027622 +0.034562 -0.005390 + H4 -0.006217 +0.004919 -0.020537 + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-FORCE (eV/Angstrom) + + ><><><><><><><><><><><><><><><><><><><><><>< + + atom x y z + C1 -0.154995 -1.114764 +0.083421 + H1 +0.292544 +0.119474 +0.337774 + H2 -0.687903 -0.019808 +0.245423 + H3 +0.710191 +0.888630 -0.138594 + H4 -0.159837 +0.126468 -0.528024 + + + ><><><><><><><><><><><><><><><><><><><><><>< + + TOTAL-STRESS (KBAR) + + ><><><><><><><><><><><><><><><><><><><><><>< + + +17.95102076 +5.24181029 -4.35918400 + +5.24181029 +13.26034469 +0.17341933 + -4.35918400 +0.17341933 -1.90016343 + + + -------------------------------------------- + !FINAL_ETOT_IS -219.6499140427659142 eV + -------------------------------------------- + + + + + + + |CLASS_NAME---------|NAME---------------|TIME(Sec)-----|CALLS----|AVG------|PER%------- + A DC_Driv reading +0.104 1 +0.10 +0.59% + A DC_Driv divide_frag +0.19 1 +0.19 +1.09% + B PW_Basis gen_pw +0.19 1 +0.19 +1.09% + A DC_Driv solve_eachf +17.20 1 +17.20 +98.32% + B Run_Frag frag_pw_line +17.20 1 +17.20 +98.32% + X FFT FFT3D +9.69 1332 +0.01 +55.38% + E potential v_of_rho +3.20 14 +0.23 +18.27% + C wavefunc wfcinit +0.22 1 +0.22 +1.24% + G Hamilt_PW cinitcgg +2.31 14 +0.17 +13.20% + H Hamilt_PW h_psi +9.29 513 +0.02 +53.10% + I Hamilt_PW add_vuspsi +0.45 513 +0.00 +2.57% + C Ions opt_ions_pw +16.63 1 +16.63 +95.07% + D electrons self_consistent +14.90 1 +14.90 +85.15% + E electrons c_bands +10.25 13 +0.79 +58.60% + F Hamilt diago +10.08 13 +0.78 +57.59% + G Diago_CG diag +7.93 13 +0.61 +45.35% + E Charge mix_rho +0.35 13 +0.03 +1.98% + ---------------------------------------------------------------------------------------- + + CLASS_NAME---------|NAME---------------|MEMORY(MB)-------- + +29.4953 + PW_Basis struc_fac +4.1212 + Use_FFT porter +4.0000 + wavefunc evc +2.0604 + Charge rho +2.0000 + Charge rho_save +2.0000 + Charge rho_core +2.0000 + potential vltot +2.0000 + potential vr +2.0000 + potential vrs +2.0000 + potential vrs1 +2.0000 + potential vnew +2.0000 + Charge rhog +1.0303 + Charge rhog_save +1.0303 + Charge rhog_core +1.0303 + ---------------------------------------------------------- + + Start Time : Fri May 7 16:18:50 2021 + Finish Time : Fri May 7 16:19:08 2021 + Total Time : +0 h +0 mins +18 secs diff --git a/tests/poscars/poscar_ref_oh.py b/tests/poscars/poscar_ref_oh.py index 23bab644f..9b12c1511 100644 --- a/tests/poscars/poscar_ref_oh.py +++ b/tests/poscars/poscar_ref_oh.py @@ -27,8 +27,12 @@ def test_cell(self): places = 6, msg = 'cell[%d][%d] failed' % (ii,jj)) - def test_frame(self): - ovito_posis = np.array([[0, 0, 0], + def test_frame(self): + if hasattr(self, "unwrap") and self.unwrap is True: + ovito_posis = np.array([[5.0739861, 2.7916155, 2.2254033], + [6.3361717, 3.4934183, 2.7767918]]) + else: + ovito_posis = np.array([[0, 0, 0], [1.2621856, 0.7018028, 0.5513885]]) for ii in range(2) : for jj in range(3) : diff --git a/tests/test_abacus_pw_scf.py b/tests/test_abacus_pw_scf.py index 3a2da17f1..f2a2c89ba 100644 --- a/tests/test_abacus_pw_scf.py +++ b/tests/test_abacus_pw_scf.py @@ -1,6 +1,6 @@ import os import numpy as np -import unittest +import unittest,shutil from context import dpdata from dpdata.unit import LengthConversion @@ -112,9 +112,30 @@ def test_energy(self) : class TestABACUSLabeledOutput(unittest.TestCase, TestABACUSSinglePointEnergy): def setUp(self): + shutil.copy('abacus.scf/INPUT.ok','abacus.scf/INPUT') self.system_ch4 = dpdata.LabeledSystem('abacus.scf',fmt='abacus/scf') # self.system_h2o = dpdata.LabeledSystem('qe.scf/02.out',fmt='qe/pw/scf') self.system_ch4_unlabeled = dpdata.System('abacus.scf/STRU.ch4', fmt='abacus/stru') + def tearDown(self): + if os.path.isfile("abacus.scf/INPUT"): + os.remove("abacus.scf/INPUT") + + +class TestABACUSLabeledOutputFail(unittest.TestCase): + + def setUp(self): + shutil.copy('abacus.scf/INPUT.fail','abacus.scf/INPUT') + self.system_ch4 = dpdata.LabeledSystem('abacus.scf',fmt='abacus/scf') + # self.system_h2o = dpdata.LabeledSystem('qe.scf/02.out',fmt='qe/pw/scf') + self.system_ch4_unlabeled = dpdata.System('abacus.scf/STRU.ch4', fmt='abacus/stru') + def tearDown(self): + if os.path.isfile("abacus.scf/INPUT"): + os.remove("abacus.scf/INPUT") + def test_return_zero(self): + self.assertEqual(len(self.system_ch4),0) + + + if __name__ == '__main__': unittest.main() diff --git a/tests/test_lammps_dump_unfold.py b/tests/test_lammps_dump_unfold.py index 4ffd5f730..68aa2c8b7 100644 --- a/tests/test_lammps_dump_unfold.py +++ b/tests/test_lammps_dump_unfold.py @@ -21,6 +21,16 @@ def test_nframes (self) : self.assertEqual(self.tmp_system.get_nframes(), 2) +class TestDumpUnwrap(unittest.TestCase, TestPOSCARoh): + def setUp(self): + self.unwrap = True + self.system = dpdata.System( + os.path.join('poscars', 'conf_unfold.dump'), + type_map=['O', 'H'], + unwrap=self.unwrap, + ) + + if __name__ == '__main__': unittest.main() diff --git a/tests/test_predict.py b/tests/test_predict.py index 8535d02ea..3ba62ec23 100644 --- a/tests/test_predict.py +++ b/tests/test_predict.py @@ -78,7 +78,7 @@ def setUp(self) : @unittest.skipIf(skip_ase,"skip ase related test. install ase to fix") -class TestASEtraj1(unittest.TestCase, CompLabeledSys, IsPBC): +class TestASEDriver(unittest.TestCase, CompLabeledSys, IsPBC): def setUp (self) : ori_sys = dpdata.LabeledSystem('poscars/deepmd.h2o.md', fmt = 'deepmd/raw', @@ -90,3 +90,34 @@ def setUp (self) : self.e_places = 6 self.f_places = 6 self.v_places = 4 + + +@unittest.skipIf(skip_ase, "skip ase related test. install ase to fix") +class TestMinimize(unittest.TestCase, CompLabeledSys, IsPBC): + def setUp (self) : + ori_sys = dpdata.LabeledSystem('poscars/deepmd.h2o.md', + fmt = 'deepmd/raw', + type_map = ['O', 'H']) + zero_driver = ZeroDriver() + self.system_1 = ori_sys.predict(driver=zero_driver) + self.system_2 = ori_sys.minimize(driver=zero_driver, minimizer="ase") + self.places = 6 + self.e_places = 6 + self.f_places = 6 + self.v_places = 4 + + +@unittest.skipIf(skip_ase, "skip ase related test. install ase to fix") +class TestMinimizeMultiSystems(unittest.TestCase, CompLabeledSys, IsPBC): + def setUp (self) : + ori_sys = dpdata.LabeledSystem('poscars/deepmd.h2o.md', + fmt = 'deepmd/raw', + type_map = ['O', 'H']) + multi_sys = dpdata.MultiSystems(ori_sys) + zero_driver = ZeroDriver() + self.system_1 = list(multi_sys.predict(driver=zero_driver).systems.values())[0] + self.system_2 = list(multi_sys.minimize(driver=zero_driver, minimizer="ase").systems.values())[0] + self.places = 6 + self.e_places = 6 + self.f_places = 6 + self.v_places = 4