|
| 1 | +import os,sys |
| 2 | +import numpy as np |
| 3 | +from .scf import ry2ev, kbar2evperang3, get_block, get_geometry_in, get_cell, get_coords |
| 4 | + |
| 5 | +# Read in geometries from an ABACUS MD trajectory. |
| 6 | +# The atomic coordinates are read in from generated files in OUT.XXXX. |
| 7 | +# Energies, forces |
| 8 | +# IMPORTANT: the program defaultly takes STRU input file as standard cell information, |
| 9 | +# therefore the direct and cartesan coordinates read could be different from the ones in |
| 10 | +# the output cif files!!! |
| 11 | +# It is highly recommanded to use ORTHOGANAL coordinates in STRU file if you wish to get |
| 12 | +# same coordinates in both dpdata and output cif files. |
| 13 | + |
| 14 | +def get_path_out(fname, inlines): |
| 15 | + # This function is different from the same-name function in scf.py. |
| 16 | + # This function returns OUT.XXXX's base directory. |
| 17 | + path_out = os.path.join(fname, "OUT.ABACUS/") |
| 18 | + for line in inlines: |
| 19 | + if len(line)>0 and "suffix" in line and "suffix"==line.split()[0]: |
| 20 | + suffix = line.split()[1] |
| 21 | + path_out = os.path.join(fname, "OUT.%s/" % suffix) |
| 22 | + break |
| 23 | + return path_out |
| 24 | + |
| 25 | +def get_coord_dump_freq(inlines): |
| 26 | + for line in inlines: |
| 27 | + if len(line)>0 and "md_dumpmdfred" in line and "md_dumpmdfred" == line.split()[0]: |
| 28 | + return int(line.split()[1]) |
| 29 | + return 1 |
| 30 | + |
| 31 | +# set up a cell according to cell info in cif file. |
| 32 | +# maybe useful later |
| 33 | +''' |
| 34 | +def setup_cell(a, b, c, alpha, beta, gamma): |
| 35 | + cell = np.zeros(3, 3) |
| 36 | + cell[0, 0] = a |
| 37 | + cell[1, 0] = b*np.cos(gamma/180*np.pi) |
| 38 | + cell[1, 1] = b*np.sin(gamma/180*np.pi) |
| 39 | + cell[2, 0] = c*np.cos(beta/180*np.pi) |
| 40 | + cell[2, 1] = c*(b*np.cos(alpha/180*np.pi) - cell[1, 0]*np.cos(beta/180*np.pi))/cell[1, 1] |
| 41 | + cell[2, 2] = np.sqrt(c**2 - cell[2, 0]**2 - cell[2, 1]**2) |
| 42 | + return cell |
| 43 | +''' |
| 44 | + |
| 45 | +def get_single_coord_from_cif(pos_file, atom_names, natoms, cell): |
| 46 | + assert(len(atom_names) == len(natoms)) |
| 47 | + nele = len(atom_names) |
| 48 | + total_natoms = sum(natoms) |
| 49 | + coord = np.zeros([total_natoms, 3]) |
| 50 | + a = 0 |
| 51 | + b = 0 |
| 52 | + c = 0 |
| 53 | + alpha = 0 |
| 54 | + beta = 0 |
| 55 | + gamma = 0 |
| 56 | + with open(pos_file, "r") as fp: |
| 57 | + lines = fp.read().split("\n") |
| 58 | + for line in lines: |
| 59 | + if "_cell_length_a" in line: |
| 60 | + a = float(line.split()[1]) |
| 61 | + if "_cell_length_b" in line: |
| 62 | + b = float(line.split()[1]) |
| 63 | + if "_cell_length_c" in line: |
| 64 | + c = float(line.split()[1]) |
| 65 | + if "_cell_angle_alpha" in line: |
| 66 | + alpha = float(line.split()[1]) |
| 67 | + if "_cell_angle_beta" in line: |
| 68 | + beta = float(line.split()[1]) |
| 69 | + if "_cell_angle_gamma" in line: |
| 70 | + gamma = float(line.split()[1]) |
| 71 | + assert(a > 0 and b > 0 and c > 0 and alpha > 0 and beta > 0 and gamma > 0) |
| 72 | + #cell = setup_cell(a, b, c, alpha, beta, gamma) |
| 73 | + coord_lines = get_block(lines=lines, keyword="_atom_site_fract_z", skip=0, nlines = total_natoms) |
| 74 | + |
| 75 | + ia_idx = 0 |
| 76 | + for it in range(nele): |
| 77 | + for ia in range(natoms[it]): |
| 78 | + coord_line = coord_lines[ia_idx].split() |
| 79 | + assert(coord_line[0] == atom_names[it]) |
| 80 | + coord[ia_idx, 0] = float(coord_line[1]) |
| 81 | + coord[ia_idx, 1] = float(coord_line[2]) |
| 82 | + coord[ia_idx, 2] = float(coord_line[3]) |
| 83 | + ia_idx+=1 |
| 84 | + coord = np.matmul(coord, cell) |
| 85 | + # important! Coordinates are converted to Cartesian coordinate. |
| 86 | + return coord |
| 87 | + |
| 88 | + |
| 89 | +def get_coords_from_cif(ndump, dump_freq, atom_names, natoms, types, path_out, cell): |
| 90 | + total_natoms = sum(natoms) |
| 91 | + #cell = np.zeros(ndump, 3, 3) |
| 92 | + coords = np.zeros([ndump, total_natoms, 3]) |
| 93 | + pos_file = os.path.join(path_out, "STRU_READIN_ADJUST.cif") |
| 94 | + # frame 0 file is different from any other frames |
| 95 | + coords[0] = get_single_coord_from_cif(pos_file, atom_names, natoms, cell) |
| 96 | + for dump_idx in range(1, ndump): |
| 97 | + pos_file = os.path.join(path_out, "md_pos_%d.cif" %(dump_idx*dump_freq)) |
| 98 | + #print("dump_idx = %s" %dump_idx) |
| 99 | + coords[dump_idx] = get_single_coord_from_cif(pos_file, atom_names, natoms, cell) |
| 100 | + return coords |
| 101 | + |
| 102 | +def get_energy_force_stress(outlines, inlines, dump_freq, ndump, natoms, atom_names): |
| 103 | + stress = None |
| 104 | + total_natoms = sum(natoms) |
| 105 | + for line in inlines: |
| 106 | + if len(line)>0 and "stress" in line and "stress" == line.split()[0] and "1" == line.split()[1]: |
| 107 | + stress = np.zeros([ndump, 3, 3]) |
| 108 | + break |
| 109 | + if type(stress) != np.ndarray: |
| 110 | + print("The ABACUS program has no stress output. Stress will not be read.") |
| 111 | + nenergy = 0 |
| 112 | + nforce = 0 |
| 113 | + nstress = 0 |
| 114 | + energy = np.zeros(ndump) |
| 115 | + force = np.zeros([ndump, total_natoms, 3]) |
| 116 | + |
| 117 | + for line_idx, line in enumerate(outlines): |
| 118 | + if "final etot is" in line: |
| 119 | + if nenergy%dump_freq == 0: |
| 120 | + energy[int(nenergy/dump_freq)] = float(line.split()[-2]) |
| 121 | + nenergy+=1 |
| 122 | + if "TOTAL-FORCE (eV/Angstrom)" in line: |
| 123 | + for iatom in range(0, total_natoms): |
| 124 | + force_line = outlines[line_idx+5+iatom] |
| 125 | + atom_force = [float(i) for i in force_line.split()[1:]] |
| 126 | + assert(len(atom_force) == 3) |
| 127 | + atom_force = np.array(atom_force) |
| 128 | + if nforce%dump_freq == 0: |
| 129 | + force[int(nforce/dump_freq), iatom] = atom_force |
| 130 | + nforce+=1 |
| 131 | + assert(nforce==nenergy) |
| 132 | + if "TOTAL-STRESS (KBAR)" in line: |
| 133 | + for idx in range(0, 3): |
| 134 | + stress_line = outlines[line_idx+4+idx] |
| 135 | + single_stress = [float(i) for i in stress_line.split()] |
| 136 | + if len(single_stress) != 3: |
| 137 | + print(single_stress) |
| 138 | + assert(len(single_stress) == 3) |
| 139 | + single_stress = np.array(single_stress) |
| 140 | + if nstress%dump_freq == 0: |
| 141 | + stress[int(nstress/dump_freq), idx] = single_stress |
| 142 | + nstress+=1 |
| 143 | + assert(nstress==nforce) |
| 144 | + if type(stress) == np.ndarray: |
| 145 | + stress *= kbar2evperang3 |
| 146 | + return energy, force, stress |
| 147 | + |
| 148 | + |
| 149 | +def get_frame (fname): |
| 150 | + if type(fname) == str: |
| 151 | + # if the input parameter is only one string, it is assumed that it is the |
| 152 | + # base directory containing INPUT file; |
| 153 | + path_in = os.path.join(fname, "INPUT") |
| 154 | + else: |
| 155 | + raise RuntimeError('invalid input') |
| 156 | + with open(path_in, 'r') as fp: |
| 157 | + inlines = fp.read().split('\n') |
| 158 | + geometry_path_in = get_geometry_in(fname, inlines) # base dir of STRU |
| 159 | + path_out = get_path_out(fname, inlines) |
| 160 | + |
| 161 | + with open(geometry_path_in, 'r') as fp: |
| 162 | + geometry_inlines = fp.read().split('\n') |
| 163 | + celldm, cell = get_cell(geometry_inlines) |
| 164 | + atom_names, natoms, types, coords = get_coords(celldm, cell, geometry_inlines, inlines) |
| 165 | + # This coords is not to be used. |
| 166 | + dump_freq = get_coord_dump_freq(inlines = inlines) |
| 167 | + ndump = int(os.popen("ls --l %s | grep 'md_pos_' | wc -l" %path_out).readlines()[0]) |
| 168 | + # number of dumped geometry files |
| 169 | + coords = get_coords_from_cif(ndump, dump_freq, atom_names, natoms, types, path_out, cell) |
| 170 | + |
| 171 | + # TODO: Read in energies, forces and pressures. |
| 172 | + with open(os.path.join(path_out, "running_md.log"), 'r') as fp: |
| 173 | + outlines = fp.read().split('\n') |
| 174 | + energy, force, stress = get_energy_force_stress(outlines, inlines, dump_freq, ndump, natoms, atom_names) |
| 175 | + if type(stress) == np.ndarray: |
| 176 | + stress *= np.linalg.det(cell) |
| 177 | + data = {} |
| 178 | + data['atom_names'] = atom_names |
| 179 | + data['atom_numbs'] = natoms |
| 180 | + data['atom_types'] = types |
| 181 | + data['cells'] = np.zeros([ndump, 3, 3]) |
| 182 | + for idx in range(ndump): |
| 183 | + data['cells'][:, :, :] = cell |
| 184 | + data['coords'] = coords |
| 185 | + data['energies'] = energy |
| 186 | + data['forces'] = force |
| 187 | + data['virials'] = stress |
| 188 | + data['orig'] = np.zeros(3) |
| 189 | + |
| 190 | + return data |
0 commit comments