11import os
2- import re
3- import sys
42import warnings
5- from ast import dump
63
74import numpy as np
85
96from .scf import (
107 bohr2ang ,
11- get_block ,
128 get_cell ,
139 get_coords ,
1410 get_geometry_in ,
1511 kbar2evperang3 ,
16- ry2ev ,
1712)
1813
1914# Read in geometries from an ABACUS MD trajectory.
@@ -72,69 +67,51 @@ def get_coords_from_dump(dumplines, natoms):
7267 for iline in range (nlines ):
7368 if "MDSTEP" in dumplines [iline ]:
7469 # read in LATTICE_CONSTANT
75- celldm = float (dumplines [iline + 1 ].split (" " )[- 1 ])
70+ # for abacus version >= v3.1.4, the unit is angstrom, and "ANGSTROM" is added at the end
71+ # for abacus version < v3.1.4, the unit is bohr
72+ celldm = float (dumplines [iline + 1 ].split ()[1 ])
73+ newversion = True
74+ if "Angstrom" not in dumplines [iline + 1 ]:
75+ celldm *= bohr2ang # transfer unit to ANGSTROM
76+ newversion = False
77+
7678 # read in LATTICE_VECTORS
7779 for ix in range (3 ):
7880 cells [iframe , ix ] = (
79- np .array (
80- [
81- float (i )
82- for i in re .split ("\s+" , dumplines [iline + 3 + ix ])[- 3 :]
83- ]
84- )
81+ np .array ([float (i ) for i in dumplines [iline + 3 + ix ].split ()[0 :3 ]])
8582 * celldm
8683 )
8784 if calc_stress :
8885 stresses [iframe , ix ] = np .array (
89- [
90- float (i )
91- for i in re .split ("\s+" , dumplines [iline + 7 + ix ])[- 3 :]
92- ]
86+ [float (i ) for i in dumplines [iline + 7 + ix ].split ()[0 :3 ]]
9387 )
88+
89+ if calc_stress :
90+ skipline = 11
91+ else :
92+ skipline = 7
93+
9494 for iat in range (total_natoms ):
95- if calc_stress :
96- coords [iframe , iat ] = (
97- np .array (
98- [
99- float (i )
100- for i in re .split ("\s+" , dumplines [iline + 11 + iat ])[
101- - 6 :- 3
102- ]
103- ]
104- )
105- * celldm
106- )
107- forces [iframe , iat ] = np .array (
108- [
109- float (i )
110- for i in re .split ("\s+" , dumplines [iline + 11 + iat ])[- 3 :]
111- ]
112- )
113- else :
114- coords [iframe , iat ] = (
115- np .array (
116- [
117- float (i )
118- for i in re .split ("\s+" , dumplines [iline + 7 + iat ])[
119- - 6 :- 3
120- ]
121- ]
122- )
123- * celldm
124- )
125- forces [iframe , iat ] = np .array (
126- [
127- float (i )
128- for i in re .split ("\s+" , dumplines [iline + 7 + iat ])[- 3 :]
129- ]
130- )
95+ # INDEX LABEL POSITION (Angstrom) FORCE (eV/Angstrom) VELOCITY (Angstrom/fs)
96+ # 0 Sn 0.000000000000 0.000000000000 0.000000000000 -0.000000000000 -0.000000000001 -0.000000000001 0.001244557166 -0.000346684288 0.000768457739
97+ # 1 Sn 0.000000000000 3.102800034079 3.102800034079 -0.000186795145 -0.000453823768 -0.000453823768 0.000550996187 -0.000886442775 0.001579501983
98+ # for abacus version >= v3.1.4, the value of POSITION is the real cartessian position, and unit is angstrom, and if cal_force the VELOCITY is added at the end.
99+ # for abacus version < v3.1.4, the real position = POSITION * celldm
100+ coords [iframe , iat ] = np .array (
101+ [float (i ) for i in dumplines [iline + skipline + iat ].split ()[2 :5 ]]
102+ )
103+
104+ if not newversion :
105+ coords [iframe , iat ] *= celldm
106+
107+ forces [iframe , iat ] = np .array (
108+ [float (i ) for i in dumplines [iline + skipline + iat ].split ()[5 :8 ]]
109+ )
131110 iframe += 1
132111 assert iframe == nframes_dump , (
133112 "iframe=%d, nframe_dump=%d. Number of frames does not match number of lines in MD_dump."
134113 % (iframe , nframes_dump )
135114 )
136- cells *= bohr2ang
137- coords *= bohr2ang
138115 stresses *= kbar2evperang3
139116 return coords , cells , forces , stresses
140117
@@ -166,12 +143,12 @@ def get_frame(fname):
166143 path_in = os .path .join (fname , "INPUT" )
167144 else :
168145 raise RuntimeError ("invalid input" )
169- with open (path_in , "r" ) as fp :
146+ with open (path_in ) as fp :
170147 inlines = fp .read ().split ("\n " )
171148 geometry_path_in = get_geometry_in (fname , inlines ) # base dir of STRU
172149 path_out = get_path_out (fname , inlines )
173150
174- with open (geometry_path_in , "r" ) as fp :
151+ with open (geometry_path_in ) as fp :
175152 geometry_inlines = fp .read ().split ("\n " )
176153 celldm , cell = get_cell (geometry_inlines )
177154 atom_names , natoms , types , coords = get_coords (
@@ -182,11 +159,11 @@ def get_frame(fname):
182159 # ndump = int(os.popen("ls -l %s | grep 'md_pos_' | wc -l" %path_out).readlines()[0])
183160 # number of dumped geometry files
184161 # coords = get_coords_from_cif(ndump, dump_freq, atom_names, natoms, types, path_out, cell)
185- with open (os .path .join (path_out , "MD_dump" ), "r" ) as fp :
162+ with open (os .path .join (path_out , "MD_dump" )) as fp :
186163 dumplines = fp .read ().split ("\n " )
187164 coords , cells , force , stress = get_coords_from_dump (dumplines , natoms )
188165 ndump = np .shape (coords )[0 ]
189- with open (os .path .join (path_out , "running_md.log" ), "r" ) as fp :
166+ with open (os .path .join (path_out , "running_md.log" )) as fp :
190167 outlines = fp .read ().split ("\n " )
191168 energy = get_energy (outlines , ndump , dump_freq )
192169
@@ -201,7 +178,7 @@ def get_frame(fname):
201178 unconv_stru += "%d " % i
202179 ndump = len (energy )
203180 if unconv_stru != "" :
204- warnings .warn (f "Structure %s are unconverged and not collected!" % unconv_stru )
181+ warnings .warn ("Structure %s are unconverged and not collected!" % unconv_stru )
205182
206183 for iframe in range (ndump ):
207184 stress [iframe ] *= np .linalg .det (cells [iframe , :, :].reshape ([3 , 3 ]))
0 commit comments