Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 73 additions & 18 deletions dpdata/cp2k/output.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
import numpy as np
import re
from collections import OrderedDict

from scipy.constants.constants import R
from .cell import cell_to_low_triangle
from ..unit import EnergyConversion, LengthConversion, ForceConversion, PressureConversion

Expand Down Expand Up @@ -29,6 +31,8 @@ def __init__(self, log_file_name, xyz_file_name, restart=False):
self.xyz_block_generator = self.get_xyz_block_generator()
self.restart_flag = restart
self.cell=None
self.print_level=None
self.atomic_kinds = None

if self.restart_flag:
self.handle_single_log_frame(next(self.log_block_generator))
Expand All @@ -43,13 +47,14 @@ def __iter__(self):
def __next__(self):
info_dict = {}
log_info_dict = self.handle_single_log_frame(next(self.log_block_generator))
#print(log_info_dict)
xyz_info_dict = self.handle_single_xyz_frame(next(self.xyz_block_generator))
eq1 = [v1==v2 for v1,v2 in zip(log_info_dict['atom_numbs'], xyz_info_dict['atom_numbs'])]
eq2 = [v1==v2 for v1,v2 in zip(log_info_dict['atom_names'], xyz_info_dict['atom_names'])]
eq3 = [v1==v2 for v1,v2 in zip(log_info_dict['atom_types'], xyz_info_dict['atom_types'])]
assert all(eq1), (log_info_dict,xyz_info_dict,'There may be errors in the file. If it is a restart task; use restart=True')
assert all(eq2), (log_info_dict,xyz_info_dict,'There may be errors in the file. If it is a restart task; use restart=True')
assert all(eq3), (log_info_dict,xyz_info_dict,'There may be errors in the file. If it is a restart task; use restart=True')
#eq1 = [v1==v2 for v1,v2 in zip(log_info_dict['atom_numbs'], xyz_info_dict['atom_numbs'])]
#eq2 = [v1==v2 for v1,v2 in zip(log_info_dict['atom_names'], xyz_info_dict['atom_names'])]
#eq3 = [v1==v2 for v1,v2 in zip(log_info_dict['atom_types'], xyz_info_dict['atom_types'])]
#assert all(eq1), (log_info_dict,xyz_info_dict,'There may be errors in the file. If it is a restart task; use restart=True')
#assert all(eq2), (log_info_dict,xyz_info_dict,'There may be errors in the file. If it is a restart task; use restart=True')
#assert all(eq3), (log_info_dict,xyz_info_dict,'There may be errors in the file. If it is a restart task; use restart=True')
assert log_info_dict['energies']==xyz_info_dict['energies'], (log_info_dict['energies'], xyz_info_dict['energies'],'There may be errors in the file')
info_dict.update(log_info_dict)
info_dict.update(xyz_info_dict)
Expand Down Expand Up @@ -103,11 +108,18 @@ def handle_single_log_frame(self, lines):
cell_angle_pattern = re.compile(r' INITIAL CELL ANGLS\[deg\]\s+=\s+(?P<alpha>\S+)\s+(?P<beta>\S+)\s+(?P<gamma>\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<ax>\S+)\s+(?P<ay>\S+)\s+(?P<az>\S+)')
cell_b_pattern = re.compile(r' CELL\| Vector b \[angstrom\]:\s+(?P<bx>\S+)\s+(?P<by>\S+)\s+(?P<bz>\S+)')
cell_c_pattern = re.compile(r' CELL\| Vector c \[angstrom\]:\s+(?P<cx>\S+)\s+(?P<cy>\S+)\s+(?P<cz>\S+)')
force_start_pattern = re.compile(r' ATOMIC FORCES in')
force_flag=False
force_end_pattern = re.compile(r' SUM OF ATOMIC FORCES')
force_lines= []
cell_flag=0
print_level_pattern = re.compile(r' GLOBAL\| Global print level\s+(?P<print_level>\S+)')
print_level_flag = 0
atomic_kinds_pattern = re.compile(r'\s+\d+\. Atomic kind:\s+(?P<akind>\S+)')
atomic_kinds = []
for line in lines:
if force_start_pattern.match(line):
force_flag=True
Expand All @@ -131,9 +143,47 @@ def handle_single_log_frame(self, lines):
cell_beta = np.deg2rad(float(cell_angle_pattern.match(line).groupdict()['beta']))
cell_gamma = np.deg2rad(float(cell_angle_pattern.match(line).groupdict()['gamma']))
cell_flag+=1
if print_level_pattern.match(line):
print_level = print_level_pattern.match(line).groupdict()['print_level']
print_level_flag += 1
if cell_a_pattern.match(line):
cell_ax = float(cell_a_pattern.match(line).groupdict()['ax'])
cell_ay = float(cell_a_pattern.match(line).groupdict()['ay'])
cell_az = float(cell_a_pattern.match(line).groupdict()['az'])
cell_flag+=1
if cell_b_pattern.match(line):
cell_bx = float(cell_b_pattern.match(line).groupdict()['bx'])
cell_by = float(cell_b_pattern.match(line).groupdict()['by'])
cell_bz = float(cell_b_pattern.match(line).groupdict()['bz'])
cell_flag+=1
if cell_c_pattern.match(line):
cell_cx = float(cell_c_pattern.match(line).groupdict()['cx'])
cell_cy = float(cell_c_pattern.match(line).groupdict()['cy'])
cell_cz = float(cell_c_pattern.match(line).groupdict()['cz'])
cell_flag+=1

if atomic_kinds_pattern.match(line):
akind = atomic_kinds_pattern.match(line).groupdict()['akind']
atomic_kinds.append(akind)
if print_level_flag == 1:
self.print_level = print_level
if print_level == 'LOW':
raise RuntimeError("please provide cp2k output with higher print level(at least MEDIUM)")


if cell_flag == 2:
self.cell = cell_to_low_triangle(cell_A,cell_B,cell_C,
cell_alpha,cell_beta,cell_gamma)
elif cell_flag == 5:
self.cell = np.asarray(
[
[cell_ax, cell_ay, cell_az],
[cell_bx, cell_by, cell_bz],
[cell_cx, cell_cy, cell_cz]]
).astype('float32')
if atomic_kinds:
self.atomic_kinds = atomic_kinds
#print(self.atomic_kinds)
# lx = cell_A
# xy = cell_B * np.cos(cell_gamma)
# xz = cell_C * np.cos(cell_beta)
Expand All @@ -146,27 +196,32 @@ def handle_single_log_frame(self, lines):

element_index = -1
element_dict = OrderedDict()
atom_types_list = []
atom_types_idx_list = []
forces_list = []
for line in force_lines[3:]:
line_list = line.split()
if element_dict.get(line_list[2]):
element_dict[line_list[2]][1]+=1
#print(line_list)
if element_dict.get(line_list[1]):
element_dict[line_list[1]][1]+=1
else:
element_index +=1
element_dict[line_list[2]]=[element_index,1]
atom_types_list.append(element_dict[line_list[2]][0])
element_dict[line_list[1]]=[element_index,1]
atom_types_idx_list.append(element_dict[line_list[1]][0])
forces_list.append([float(line_list[3])*AU_TO_EV_EVERY_ANG,
float(line_list[4])*AU_TO_EV_EVERY_ANG,
float(line_list[5])*AU_TO_EV_EVERY_ANG])

atom_names=list(element_dict.keys())
#print(atom_types_idx_list)
#atom_names=list(element_dict.keys())
atom_names=self.atomic_kinds
atom_numbs=[]
for ii in atom_names:

for ii in element_dict.keys():
atom_numbs.append(element_dict[ii][1])
#print(atom_numbs)
info_dict['atom_names'] = atom_names
info_dict['atom_numbs'] = atom_numbs
info_dict['atom_types'] = np.asarray(atom_types_list)
info_dict['atom_types'] = np.asarray(atom_types_idx_list)
info_dict['print_level'] = self.print_level
info_dict['cells'] = np.asarray([self.cell]).astype('float32')
info_dict['energies'] = np.asarray([energy]).astype('float32')
info_dict['forces'] = np.asarray([forces_list]).astype('float32')
Expand Down Expand Up @@ -208,9 +263,9 @@ def handle_single_xyz_frame(self, lines):
atom_numbs=[]
for ii in atom_names:
atom_numbs.append(element_dict[ii][1])
info_dict['atom_names'] = atom_names
info_dict['atom_numbs'] = atom_numbs
info_dict['atom_types'] = np.asarray(atom_types_list)
#info_dict['atom_names'] = atom_names
#info_dict['atom_numbs'] = atom_numbs
#info_dict['atom_types'] = np.asarray(atom_types_list)
info_dict['coords'] = np.asarray([coords_list]).astype('float32')
info_dict['energies'] = np.array([energy]).astype('float32')
info_dict['orig']=[0,0,0]
Expand Down
31 changes: 31 additions & 0 deletions tests/cp2k/aimd/DPGEN-1.ener
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
# Step Nr. Time[fs] Kin.[a.u.] Temp[K] Pot.[a.u.] Cons Qty[a.u.] UsedTime[s]
1 0.000000 0.000000000 0.000000000 -12660.628917195 0.000000000 936.038164850
2 0.000000 0.000000000 0.000000000 -12660.653046555 0.000000000 274.698523818
3 0.000000 0.000000000 0.000000000 -12660.715916301 0.000000000 344.194269081
4 0.000000 0.000000000 0.000000000 -12660.683902467 0.000000000 365.777609413
5 0.000000 0.000000000 0.000000000 -12660.717204918 0.000000000 355.401395731
6 0.000000 0.000000000 0.000000000 -12660.733526044 0.000000000 441.982779772
7 0.000000 0.000000000 0.000000000 -12660.692762118 0.000000000 399.540830899
8 0.000000 0.000000000 0.000000000 -12660.693929365 0.000000000 441.216501789
9 0.000000 0.000000000 0.000000000 -12660.651822948 0.000000000 397.996012313
10 0.000000 0.000000000 0.000000000 -12660.695752782 0.000000000 387.826853781
11 0.000000 0.000000000 0.000000000 -12660.648232319 0.000000000 420.456175307
12 0.000000 0.000000000 0.000000000 -12660.607125395 0.000000000 546.109416880
13 0.000000 0.000000000 0.000000000 -12660.628884156 0.000000000 419.908203241
14 0.000000 0.000000000 0.000000000 -12660.608694793 0.000000000 511.207648837
15 0.000000 0.000000000 0.000000000 -12660.632201776 0.000000000 419.362851471
16 0.000000 0.000000000 0.000000000 -12660.622454799 0.000000000 430.431832543
17 0.000000 0.000000000 0.000000000 -12660.618601776 0.000000000 441.716595725
18 0.000000 0.000000000 0.000000000 -12660.670388010 0.000000000 365.645702597
19 0.000000 0.000000000 0.000000000 -12660.647822007 0.000000000 408.940343004
20 0.000000 0.000000000 0.000000000 -12660.656801404 0.000000000 451.677564803
21 0.000000 0.000000000 0.000000000 -12660.658641834 0.000000000 623.018366899
22 0.000000 0.000000000 0.000000000 -12660.693126694 0.000000000 462.518734109
23 0.000000 0.000000000 0.000000000 -12660.678342029 0.000000000 355.707985190
24 0.000000 0.000000000 0.000000000 -12660.652483744 0.000000000 398.295128799
25 0.000000 0.000000000 0.000000000 -12660.631736112 0.000000000 442.329493437
26 0.000000 0.000000000 0.000000000 -12660.641113462 0.000000000 397.965252467
27 0.000000 0.000000000 0.000000000 -12660.713753654 0.000000000 441.802233378
28 0.000000 0.000000000 0.000000000 -12660.724079194 0.000000000 397.543293045
29 0.000000 0.000000000 0.000000000 -12660.705785856 0.000000000 355.694321892
30 0.000000 0.000000000 0.000000000 -12660.703546464 0.000000000 420.815910676
Loading