【dpdata】How to Read Virials Data from cp.x Output Using dpdata?

[hellolori]

Dear Developers,
I need to read the energies, forces, and virials data from files generated by cp.x. According to the dpdata manual, which states:

qe/cp/traj: Quantum Espresso CP trajectory files should include: file_name+'.in', file_name+'.pos', file_name+'.evp', and file_name+'.for'

I have successfully read the energies and forces data using the following script:

#!/usr/bin/env python
import dpdata
import numpy as np
import os
import shutil
import argparse
data = dpdata.LabeledSystem(file_name="si", fmt="qe/cp/traj")
print(data["energies"])
print(data["forces"])

However, I do not know how to read the virials data, which appears to be written in the file file_name+'.str'. This file is not mentioned in the manual. Could you please advise me on how to read the virials data?

[wanghan-iapcm]

Loading virial tensor is not supported by qe/cp/traj. you are definitely welcome to contribute.

[hellolori]

Sure. I have modified traj.py to support loading the virial tensor from the prefix.str file (if present). The updated script can now read the stress tensor from prefix.str and the volume information from the prefix.evp file, and then convert the stress tensor to the virial tensor in units of eV. The modified traj.py is attached as traj.readstress.py.tgz.

traj.readstress.py.tgz

$ diff traj.readstress.py traj.ori.py
6,7d5
< import os
<
23,42d20
< gpa2evperbohr = 0.00092489622836616
<
< def load_evpfile_column(fname, col, begin=0, step=1):
<     """
<     Read specified column from .evp file (0-based index)
<     """
<     data = np.loadtxt(fname)
<     steps = []
<     for ii in data[begin::step, 0]:
<         steps.append("%d" % ii)
<     with open(fname) as fp:
<         while True:
<             line = fp.readline()
<             if not line:
<                 return None
<             if line.split()[0][0] != "#":
<                 nw = len(line.split())
<                 break
<     data = np.reshape(data, [-1, nw])
<     return data[begin::step, col], steps
228d205
<
248,266d224
<
<     ### New part: Read stress and volume, calculate virials ###
<     stress_fname = prefix + ".str"
<     evp_fname = prefix + ".evp"
<     if os.path.exists(stress_fname):
<         # 1. Read stress tensor (in GPa) for each structure
<         stress, ssteps = load_data(
<             stress_fname, 3, begin=begin, step=step, convert=1.0
<         )
<         # 2. Read volume from the 10th column (Python index 9) in .evp file
<         volumes, vsteps = load_evpfile_column(evp_fname, 9, begin=begin, step=step)
<         volumes = np.array(volumes).reshape(-1)
<         # 3. Calculate virials for each structure
<         virials = []
<         for i in range(stress.shape[0]):
<             virials.append(stress[i] * gpa2evperbohr * volumes[i])
<         data["virials"] = np.array(virials)  # shape: [nframe, 3, 3]
<     ### End of new part ###

[wanghan-iapcm]

That's great! thank you very much
could you please provide us an example set of output files and the expected stress value, so we may add a unit test.
Note that we are not requesting the whole output of production runs, but a few lines from them so we may test if the code works as expected.

[hellolori]

Sure. Please find attached the input and output files from cp.x calculations for Si. In addition to the files required by dpdata, I have included a .str file, which stores the stress tensors (in GPa) for each structure.
Using the following command:

data = dpdata.LabeledSystem(file_name="si", fmt="qe/cp/traj")

the presence of the .str file will be detected, and the stress tensors will be converted to virial tensors with units of eV, and stored as data[0:num_of_structure-1]["virials"].
Please let me know the results of your tests, so I can confirm that the modified traj.py is working as expected.

si_cpmd.tgz