Skip to content

parse virial info in cp2k #113

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
amcadmus merged 1 commit into from
Aug 14, 2020
Merged

parse virial info in cp2k #113

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
31 changes: 29 additions & 2 deletions dpdata/cp2k/output.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -210,13 +210,16 @@ def handle_single_xyz_frame(self, lines):
def get_frames (fname) :
coord_flag = False
force_flag = False
stress_flag = False
eV = 2.72113838565563E+01 # hatree to eV
angstrom = 5.29177208590000E-01 # Bohrto Angstrom
angstrom = 5.29177208590000E-01 # Bohr to Angstrom
GPa = 160.21766208 # 1 eV/(Angstrom^3) = 160.21 GPa
fp = open(fname)
atom_symbol_list = []
cell = []
coord = []
force = []
stress = []
coord_count = 0
for idx, ii in enumerate(fp) :
if 'CELL| Vector' in ii :
Expand Down Expand Up @@ -244,10 +247,23 @@ def get_frames (fname) :
force_flag = False
else :
force.append(ii.split()[3:6])
# add reading stress tensor
if 'STRESS TENSOR [GPa' in ii :
stress_flag = True
stress_idx = idx
if stress_flag :
if (idx > stress_idx + 2):
if (ii == '\n') :
stress_flag = False
else :
stress.append(ii.split()[1:4])


fp.close()
assert(coord), "cannot find coords"
assert(energy), "cannot find energies"
assert(force), "cannot find forces"
assert(stress), "cannot find stress"

#conver to float array and add extra dimension for nframes
cell = np.array(cell)
Expand All @@ -260,10 +276,21 @@ def get_frames (fname) :
force = np.array(force)
force = force.astype(np.float)
force = force[np.newaxis, :, :]
stress = np.array(stress)
stress = stress.astype(np.float)
stress = stress[np.newaxis, :, :]
# force unit conversion, default unit in cp2k is hartree/bohr
force = force * eV / angstrom
# energy unit conversion, default unit in cp2k is hartree
energy = float(energy) * eV
energy = np.array(energy)
energy = energy[np.newaxis]
# stress to virial conversion, default unit in cp2k is GPa
# note the stress is virial = stress * volume
virial = stress * np.linalg.det(cell[0])/GPa



tmp_names, symbol_idx = np.unique(atom_symbol_list, return_index=True)
atom_types = []
atom_numbs = []
Expand All @@ -278,7 +305,7 @@ def get_frames (fname) :

atom_types = np.array(atom_types)

return list(atom_names), atom_numbs, atom_types, cell, coord, energy, force
return list(atom_names), atom_numbs, atom_types, cell, coord, energy, force, virial



Expand Down
61 changes: 31 additions & 30 deletions dpdata/system.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ def from_fmt(self, file_name, fmt='auto', **kwargs):
func(self, file_name, **kwargs)
else :
raise RuntimeError('unknow data format ' + fmt)

def to(self, fmt, *args, **kwargs):
fmt = fmt.lower()
to_funcs = self.register_to_funcs.funcs
Expand Down Expand Up @@ -349,15 +349,15 @@ def append(self, system) :
for ii in ['coords', 'cells'] :
self.data[ii] = np.concatenate((self.data[ii], system[ii]), axis = 0)
if self.nopbc and not system.nopbc:
# appended system uses PBC, cancel nopbc
# appended system uses PBC, cancel nopbc
self.data['nopbc'] = False
return True

def sort_atom_names(self, type_map=None):
"""
Sort atom_names of the system and reorder atom_numbs and atom_types accoarding
to atom_names. If type_map is not given, atom_names will be sorted by
alphabetical order. If type_map is given, atom_names will be type_map.
alphabetical order. If type_map is given, atom_names will be type_map.

Parameters
----------
Expand Down Expand Up @@ -394,7 +394,7 @@ def check_type_map(self, type_map):
Parameters
----------
type_map : list
type_map
type_map
"""
if type_map is not None and type_map != self.data['atom_names']:
self.sort_atom_names(type_map=type_map)
Expand All @@ -416,13 +416,13 @@ def formula(self):
"""
Return the formula of this system, like C3H5O2
"""
return ''.join(["{}{}".format(symbol,numb) for symbol,numb in
return ''.join(["{}{}".format(symbol,numb) for symbol,numb in
zip(self.data['atom_names'], self.data['atom_numbs'])])

@property
def uniq_formula(self):
"""
Return the uniq_formula of this system.
Return the uniq_formula of this system.
The uniq_formula sort the elements in formula by names.
Systems with the same uniq_formula can be append together.
"""
Expand Down Expand Up @@ -456,9 +456,9 @@ def apply_pbc(self) :
def remove_pbc(self, protect_layer = 9):
"""
This method does NOT delete the definition of the cells, it
(1) revises the cell to a cubic cell and ensures that the cell
(1) revises the cell to a cubic cell and ensures that the cell
boundary to any atom in the system is no less than `protect_layer`
(2) translates the system such that the center-of-geometry of the system
(2) translates the system such that the center-of-geometry of the system
locates at the center of the cell.

Parameters
Expand All @@ -474,12 +474,12 @@ def remove_pbc(self, protect_layer = 9):
cog = np.average(tmpcoord, axis = 0)
dist = tmpcoord - np.tile(cog, [natoms, 1])
max_dist = np.max(np.linalg.norm(dist, axis = 1))
h_cell_size = max_dist + protect_layer
h_cell_size = max_dist + protect_layer
cell_size = h_cell_size * 2
shift = np.array([1,1,1]) * h_cell_size - cog
self.data['coords'][ff] = self.data['coords'][ff] + np.tile(shift, [natoms, 1])
self.data['cells'][ff] = cell_size * np.eye(3)
self.data['cells'][ff] = cell_size * np.eye(3)


@register_from_funcs.register_funcs("lmp")
@register_from_funcs.register_funcs("lammps/lmp")
Expand Down Expand Up @@ -508,7 +508,7 @@ def to_pymatgen_structure(self):
structure=Structure(system.data['cells'][0],species,system.data['coords'][0],coords_are_cartesian=True)
structures.append(structure)
return structures


@register_to_funcs.register_funcs("ase/structure")
def to_ase_structure(self):
Expand Down Expand Up @@ -570,7 +570,7 @@ def from_vasp_poscar(self, file_name) :
def to_vasp_string(self, frame_idx=0):
"""
Dump the system in vasp POSCAR format string

Parameters
----------
frame_idx : int
Expand All @@ -579,7 +579,7 @@ def to_vasp_string(self, frame_idx=0):
assert(frame_idx < len(self.data['coords']))
w_str = dpdata.vasp.poscar.from_system_data(self.data, frame_idx)
return w_str

@register_to_funcs.register_funcs("vasp/poscar")
def to_vasp_poscar(self, file_name, frame_idx = 0) :
"""
Expand Down Expand Up @@ -660,7 +660,7 @@ def to_deepmd_raw(self, folder) :
"""
Dump the system in deepmd raw format to `folder`
"""
dpdata.deepmd.raw.dump(folder, self.data)
dpdata.deepmd.raw.dump(folder, self.data)

@register_from_funcs.register_funcs('siesta/output')
def from_siesta_output(self, fname):
Expand Down Expand Up @@ -761,7 +761,7 @@ def replicate(self, ncopy):

Parameters
----------
ncopy :
ncopy :
list: [4,2,3]
or tuple: (4,2,3,)
make `ncopy[0]` copys in x dimensions,
Expand All @@ -778,7 +778,7 @@ def replicate(self, ncopy):
for ii in ncopy:
if type(ii) is not int:
raise RuntimeError('ncopy must be a list or tuple must with 3 int')

tmp = System()
nframes = self.get_nframes()
data = self.data
Expand All @@ -799,10 +799,10 @@ def replicate(self, ncopy):
tmp.data['coords'] = np.reshape(np.transpose(tmp.data['coords'], [3,4,0,1,2,5]), (nframes, -1 , 3))
return tmp

def perturb(self,
pert_num,
def perturb(self,
pert_num,
cell_pert_fraction,
atom_pert_distance,
atom_pert_distance,
atom_pert_style='normal'):
"""
Perturb each frame in the system randomly.
Expand All @@ -816,8 +816,8 @@ def perturb(self,
That means the system to be returned will contain `pert_num` * frame_num of the input system.
cell_pert_fraction : float
A fraction determines how much (relatively) will cell deform.
The cell of each frame is deformed by a symmetric matrix perturbed from identity.
The perturbation to the diagonal part is subject to a uniform distribution in [-cell_pert_fraction, cell_pert_fraction),
The cell of each frame is deformed by a symmetric matrix perturbed from identity.
The perturbation to the diagonal part is subject to a uniform distribution in [-cell_pert_fraction, cell_pert_fraction),
and the perturbation to the off-diagonal part is subject to a uniform distribution in [-0.5*cell_pert_fraction, 0.5*cell_pert_fraction).
atom_pert_distance: float
unit: Angstrom. A distance determines how far atoms will move.
Expand All @@ -839,7 +839,7 @@ def perturb(self,
The perturbed_system. It contains `pert_num` * frame_num of the input system frames.
"""
perturbed_system = System()
nframes = self.get_nframes()
nframes = self.get_nframes()
for ii in range(nframes):
for jj in range(pert_num):
tmp_system = self[ii].copy()
Expand Down Expand Up @@ -921,7 +921,7 @@ def get_atom_perturb_vector(atom_pert_distance, atom_pert_style='normal'):
random_vector = None
if atom_pert_distance < 0:
raise RuntimeError('atom_pert_distance can not be negative')

if atom_pert_style == 'normal':
e = np.random.randn(3)
random_vector=(atom_pert_distance/np.sqrt(3))*e
Expand Down Expand Up @@ -1173,7 +1173,7 @@ def from_qe_cp_traj(self, prefix, begin = 0, step = 1) :
)
self.data['energies'], self.data['forces'], es \
= dpdata.qe.traj.to_system_label(prefix + '.in', prefix, begin = begin, step = step)
assert(cs == es), "the step key between files are not consistent"
assert(cs == es), "the step key between files are not consistent"
self.rot_lower_triangular()

@register_from_funcs.register_funcs('qe/pw/scf')
Expand All @@ -1188,7 +1188,7 @@ def from_qe_pw_scf(self, file_name) :
self.data['virials'], \
= dpdata.qe.scf.get_frame(file_name)
self.rot_lower_triangular()

@register_from_funcs.register_funcs('siesta/output')
def from_siesta_output(self, file_name) :
self.data['atom_names'], \
Expand All @@ -1213,15 +1213,15 @@ def from_siesta_aiMD_output(self, file_name):
self.data['forces'], \
self.data['virials'] \
= dpdata.siesta.aiMD_output.get_aiMD_frame(file_name)

@register_from_funcs.register_funcs('gaussian/log')
def from_gaussian_log(self, file_name, md=False):
try:
self.data = dpdata.gaussian.log.to_system_data(file_name, md=md)
except AssertionError:
self.data['energies'], self.data['forces']= [], []
self.data['nopbc'] = True

@register_from_funcs.register_funcs('gaussian/md')
def from_gaussian_md(self, file_name):
self.from_gaussian_log(file_name, md=True)
Expand All @@ -1248,6 +1248,7 @@ def from_cp2k_output(self, file_name) :
self.data['coords'], \
self.data['energies'], \
self.data['forces'], \
self.data['virials'] \
= dpdata.cp2k.output.get_frames(file_name)
@register_from_funcs.register_funcs('movement')
@register_from_funcs.register_funcs('MOVEMENT')
Expand Down Expand Up @@ -1405,7 +1406,7 @@ def __add__(self, others) :
elif isinstance(others, list):
return self.__class__(self, *others)
raise RuntimeError("Unspported data structure")

@classmethod
def from_file(cls,file_name,fmt):
multi_systems = cls()
Expand Down Expand Up @@ -1547,7 +1548,7 @@ def check_System(data):
if len(data['coords']) > 0 :
assert( len(data['coords'][0])==len(data['atom_types'])==sum(data['atom_numbs']) )
else :
assert( len(data['atom_types'])==sum(data['atom_numbs']) )
assert( len(data['atom_types'])==sum(data['atom_numbs']) )
assert( len(data['cells']) == len(data['coords']) )
assert( len(data['atom_names'])==len(data['atom_numbs']) )

Expand Down
Loading