Coverage for /builds/ase/ase/ase/io/lammpsdata.py: 94.08%
338 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import re
4import warnings
6import numpy as np
8from ase.atoms import Atoms
9from ase.calculators.lammps import Prism, convert
10from ase.data import atomic_masses, atomic_numbers
11from ase.utils import reader, writer
14@reader
15def read_lammps_data(
16 fileobj,
17 Z_of_type: dict = None,
18 sort_by_id: bool = True,
19 read_image_flags: bool = True,
20 units: str = 'metal',
21 atom_style: str = None,
22 style: str = None,
23):
24 """Method which reads a LAMMPS data file.
26 Parameters
27 ----------
28 fileobj : file | str
29 File from which data should be read.
30 Z_of_type : dict[int, int], optional
31 Mapping from LAMMPS atom types (typically starting from 1) to atomic
32 numbers. If None, if there is the "Masses" section, atomic numbers are
33 guessed from the atomic masses. Otherwise, atomic numbers of 1 (H), 2
34 (He), etc. are assigned to atom types of 1, 2, etc. Default is None.
35 sort_by_id : bool, optional
36 Order the particles according to their id. Might be faster to set it
37 False. Default is True.
38 read_image_flags: bool, default True
39 If True, the lattice translation vectors derived from image flags are
40 added to atomic positions.
41 units : str, optional
42 `LAMMPS units <https://docs.lammps.org/units.html>`__. Default is
43 'metal'.
44 atom_style : {'atomic', 'charge', 'full'} etc., optional
45 `LAMMPS atom style <https://docs.lammps.org/atom_style.html>`__.
46 If None, `atom_style` is guessed in the following priority (1) comment
47 after `Atoms` (2) length of fields (valid only `atomic` and `full`).
48 Default is None.
49 """
50 if style is not None:
51 warnings.warn(
52 FutureWarning('"style" is deprecated; please use "atom_style".'),
53 )
54 atom_style = style
55 # begin read_lammps_data
56 file_comment = next(fileobj).rstrip()
58 # default values (https://docs.lammps.org/read_data.html)
59 # in most cases these will be updated below
60 natoms = 0
61 # N_types = 0
62 xlo, xhi = -0.5, 0.5
63 ylo, yhi = -0.5, 0.5
64 zlo, zhi = -0.5, 0.5
65 xy, xz, yz = 0.0, 0.0, 0.0
67 mass_in = {}
68 vel_in = {}
69 bonds_in = []
70 angles_in = []
71 dihedrals_in = []
73 sections = [
74 'Atoms',
75 'Velocities',
76 'Masses',
77 'Charges',
78 'Ellipsoids',
79 'Lines',
80 'Triangles',
81 'Bodies',
82 'Bonds',
83 'Angles',
84 'Dihedrals',
85 'Impropers',
86 'Impropers Pair Coeffs',
87 'PairIJ Coeffs',
88 'Pair Coeffs',
89 'Bond Coeffs',
90 'Angle Coeffs',
91 'Dihedral Coeffs',
92 'Improper Coeffs',
93 'BondBond Coeffs',
94 'BondAngle Coeffs',
95 'MiddleBondTorsion Coeffs',
96 'EndBondTorsion Coeffs',
97 'AngleTorsion Coeffs',
98 'AngleAngleTorsion Coeffs',
99 'BondBond13 Coeffs',
100 'AngleAngle Coeffs',
101 ]
102 header_fields = [
103 'atoms',
104 'bonds',
105 'angles',
106 'dihedrals',
107 'impropers',
108 'atom types',
109 'bond types',
110 'angle types',
111 'dihedral types',
112 'improper types',
113 'extra bond per atom',
114 'extra angle per atom',
115 'extra dihedral per atom',
116 'extra improper per atom',
117 'extra special per atom',
118 'ellipsoids',
119 'lines',
120 'triangles',
121 'bodies',
122 'xlo xhi',
123 'ylo yhi',
124 'zlo zhi',
125 'xy xz yz',
126 ]
127 sections_re = '(' + '|'.join(sections).replace(' ', '\\s+') + ')'
128 header_fields_re = '(' + '|'.join(header_fields).replace(' ', '\\s+') + ')'
130 section = None
131 header = True
132 for line in fileobj:
133 # get string after #; if # does not exist, return ''
134 line_comment = re.sub(r'^.*#|^.*$', '', line).strip()
135 line = re.sub('#.*', '', line).rstrip().lstrip()
136 if re.match('^\\s*$', line): # skip blank lines
137 continue
139 # check for known section names
140 match = re.match(sections_re, line)
141 if match is not None:
142 section = match.group(0).rstrip().lstrip()
143 header = False
144 if section == 'Atoms': # id *
145 # guess `atom_style` from the comment after `Atoms` if exists
146 if atom_style is None and line_comment != '':
147 atom_style = line_comment
148 mol_id_in, type_in, charge_in, pos_in, travel_in = \
149 _read_atoms_section(fileobj, natoms, atom_style)
150 continue
152 if header:
153 field = None
154 val = None
155 match = re.match('(.*)\\s+' + header_fields_re, line)
156 if match is not None:
157 field = match.group(2).lstrip().rstrip()
158 val = match.group(1).lstrip().rstrip()
159 if field is not None and val is not None:
160 if field == 'atoms':
161 natoms = int(val)
162 elif field == 'xlo xhi':
163 (xlo, xhi) = (float(x) for x in val.split())
164 elif field == 'ylo yhi':
165 (ylo, yhi) = (float(x) for x in val.split())
166 elif field == 'zlo zhi':
167 (zlo, zhi) = (float(x) for x in val.split())
168 elif field == 'xy xz yz':
169 (xy, xz, yz) = (float(x) for x in val.split())
171 if section is not None:
172 fields = line.split()
173 if section == 'Velocities': # id vx vy vz
174 atom_id = int(fields[0])
175 vel_in[atom_id] = [float(fields[_]) for _ in (1, 2, 3)]
176 elif section == 'Masses':
177 mass_in[int(fields[0])] = float(fields[1])
178 elif section == 'Bonds': # id type atom1 atom2
179 bonds_in.append([int(fields[_]) for _ in (1, 2, 3)])
180 elif section == 'Angles': # id type atom1 atom2 atom3
181 angles_in.append([int(fields[_]) for _ in (1, 2, 3, 4)])
182 elif section == 'Dihedrals': # id type atom1 atom2 atom3 atom4
183 dihedrals_in.append([int(fields[_]) for _ in (1, 2, 3, 4, 5)])
185 # set cell
186 cell = np.zeros((3, 3))
187 cell[0, 0] = xhi - xlo
188 cell[1, 1] = yhi - ylo
189 cell[2, 2] = zhi - zlo
190 cell[1, 0] = xy
191 cell[2, 0] = xz
192 cell[2, 1] = yz
194 # initialize arrays for per-atom quantities
195 positions = np.zeros((natoms, 3))
196 numbers = np.zeros(natoms, int)
197 ids = np.zeros(natoms, int)
198 types = np.zeros(natoms, int)
199 velocities = np.zeros((natoms, 3)) if len(vel_in) > 0 else None
200 masses = np.zeros(natoms) if len(mass_in) > 0 else None
201 mol_id = np.zeros(natoms, int) if len(mol_id_in) > 0 else None
202 charge = np.zeros(natoms, float) if len(charge_in) > 0 else None
203 travel = np.zeros((natoms, 3), int) if len(travel_in) > 0 else None
204 bonds = [''] * natoms if len(bonds_in) > 0 else None
205 angles = [''] * natoms if len(angles_in) > 0 else None
206 dihedrals = [''] * natoms if len(dihedrals_in) > 0 else None
208 ind_of_id = {}
209 # copy per-atom quantities from read-in values
210 for (i, atom_id) in enumerate(pos_in.keys()):
211 # by id
212 if sort_by_id:
213 ind = atom_id - 1
214 else:
215 ind = i
216 ind_of_id[atom_id] = ind
217 atom_type = type_in[atom_id]
218 positions[ind, :] = pos_in[atom_id]
219 if velocities is not None:
220 velocities[ind, :] = vel_in[atom_id]
221 if travel is not None:
222 travel[ind] = travel_in[atom_id]
223 if mol_id is not None:
224 mol_id[ind] = mol_id_in[atom_id]
225 if charge is not None:
226 charge[ind] = charge_in[atom_id]
227 ids[ind] = atom_id
228 # by type
229 types[ind] = atom_type
230 if Z_of_type is None:
231 numbers[ind] = atom_type
232 else:
233 numbers[ind] = Z_of_type[atom_type]
234 if masses is not None:
235 masses[ind] = mass_in[atom_type]
236 # convert units
237 positions = convert(positions, 'distance', units, 'ASE')
238 cell = convert(cell, 'distance', units, 'ASE')
239 if masses is not None:
240 masses = convert(masses, 'mass', units, 'ASE')
241 if velocities is not None:
242 velocities = convert(velocities, 'velocity', units, 'ASE')
244 # guess atomic numbers from atomic masses
245 # this must be after the above mass-unit conversion
246 if Z_of_type is None and masses is not None:
247 numbers = _masses2numbers(masses)
249 # create ase.Atoms
250 atoms = Atoms(
251 positions=positions,
252 numbers=numbers,
253 masses=masses,
254 cell=cell,
255 pbc=[True, True, True],
256 )
258 # add lattice translation vectors
259 if read_image_flags and travel is not None:
260 scaled_positions = atoms.get_scaled_positions(wrap=False) + travel
261 atoms.set_scaled_positions(scaled_positions)
263 # set velocities (can't do it via constructor)
264 if velocities is not None:
265 atoms.set_velocities(velocities)
267 atoms.arrays['id'] = ids
268 atoms.arrays['type'] = types
269 if mol_id is not None:
270 atoms.arrays['mol-id'] = mol_id
271 if charge is not None:
272 atoms.arrays['initial_charges'] = charge
273 atoms.arrays['mmcharges'] = charge.copy()
275 if bonds is not None:
276 for (atom_type, at1, at2) in bonds_in:
277 i_a1 = ind_of_id[at1]
278 i_a2 = ind_of_id[at2]
279 if len(bonds[i_a1]) > 0:
280 bonds[i_a1] += ','
281 bonds[i_a1] += f'{i_a2:d}({atom_type:d})'
282 for i, bond in enumerate(bonds):
283 if len(bond) == 0:
284 bonds[i] = '_'
285 atoms.arrays['bonds'] = np.array(bonds)
287 if angles is not None:
288 for (atom_type, at1, at2, at3) in angles_in:
289 i_a1 = ind_of_id[at1]
290 i_a2 = ind_of_id[at2]
291 i_a3 = ind_of_id[at3]
292 if len(angles[i_a2]) > 0:
293 angles[i_a2] += ','
294 angles[i_a2] += f'{i_a1:d}-{i_a3:d}({atom_type:d})'
295 for i, angle in enumerate(angles):
296 if len(angle) == 0:
297 angles[i] = '_'
298 atoms.arrays['angles'] = np.array(angles)
300 if dihedrals is not None:
301 for (atom_type, at1, at2, at3, at4) in dihedrals_in:
302 i_a1 = ind_of_id[at1]
303 i_a2 = ind_of_id[at2]
304 i_a3 = ind_of_id[at3]
305 i_a4 = ind_of_id[at4]
306 if len(dihedrals[i_a1]) > 0:
307 dihedrals[i_a1] += ','
308 dihedrals[i_a1] += f'{i_a2:d}-{i_a3:d}-{i_a4:d}({atom_type:d})'
309 for i, dihedral in enumerate(dihedrals):
310 if len(dihedral) == 0:
311 dihedrals[i] = '_'
312 atoms.arrays['dihedrals'] = np.array(dihedrals)
314 atoms.info['comment'] = file_comment
316 return atoms
319def _read_atoms_section(fileobj, natoms: int, style: str = None):
320 type_in = {}
321 mol_id_in = {}
322 charge_in = {}
323 pos_in = {}
324 travel_in = {}
325 next(fileobj) # skip blank line just after `Atoms`
326 for _ in range(natoms):
327 line = next(fileobj)
328 line = re.sub('#.*', '', line).rstrip().lstrip()
329 fields = line.split()
330 if style is None:
331 style = _guess_atom_style(fields)
332 atom_id = int(fields[0])
333 if style == 'full' and len(fields) in (7, 10):
334 # id mol-id type q x y z [tx ty tz]
335 type_in[atom_id] = int(fields[2])
336 pos_in[atom_id] = tuple(float(fields[_]) for _ in (4, 5, 6))
337 mol_id_in[atom_id] = int(fields[1])
338 charge_in[atom_id] = float(fields[3])
339 if len(fields) == 10:
340 travel_in[atom_id] = tuple(int(fields[_]) for _ in (7, 8, 9))
341 elif style == 'atomic' and len(fields) in (5, 8):
342 # id type x y z [tx ty tz]
343 type_in[atom_id] = int(fields[1])
344 pos_in[atom_id] = tuple(float(fields[_]) for _ in (2, 3, 4))
345 if len(fields) == 8:
346 travel_in[atom_id] = tuple(int(fields[_]) for _ in (5, 6, 7))
347 elif style in ('angle', 'bond', 'molecular') and len(fields) in (6, 9):
348 # id mol-id type x y z [tx ty tz]
349 type_in[atom_id] = int(fields[2])
350 pos_in[atom_id] = tuple(float(fields[_]) for _ in (3, 4, 5))
351 mol_id_in[atom_id] = int(fields[1])
352 if len(fields) == 9:
353 travel_in[atom_id] = tuple(int(fields[_]) for _ in (6, 7, 8))
354 elif style == 'charge' and len(fields) in (6, 9):
355 # id type q x y z [tx ty tz]
356 type_in[atom_id] = int(fields[1])
357 pos_in[atom_id] = tuple(float(fields[_]) for _ in (3, 4, 5))
358 charge_in[atom_id] = float(fields[2])
359 if len(fields) == 9:
360 travel_in[atom_id] = tuple(int(fields[_]) for _ in (6, 7, 8))
361 else:
362 raise RuntimeError(
363 f'Style "{style}" not supported or invalid. '
364 f'Number of fields: {len(fields)}'
365 )
366 return mol_id_in, type_in, charge_in, pos_in, travel_in
369def _guess_atom_style(fields):
370 """Guess `atom_sytle` from the length of fields."""
371 if len(fields) in (5, 8):
372 return 'atomic'
373 if len(fields) in (7, 10):
374 return 'full'
375 raise ValueError('atom_style cannot be guessed from len(fields)')
378def _masses2numbers(masses):
379 """Guess atomic numbers from atomic masses."""
380 return np.argmin(np.abs(atomic_masses - masses[:, None]), axis=1)
383@writer
384def write_lammps_data(
385 fd,
386 atoms: Atoms,
387 *,
388 specorder: list = None,
389 reduce_cell: bool = False,
390 force_skew: bool = False,
391 prismobj: Prism = None,
392 write_image_flags: bool = False,
393 masses: bool = False,
394 velocities: bool = False,
395 units: str = 'metal',
396 bonds: bool = True,
397 atom_style: str = 'atomic',
398):
399 """Write atomic structure data to a LAMMPS data file.
401 Parameters
402 ----------
403 fd : file|str
404 File to which the output will be written.
405 atoms : Atoms
406 Atoms to be written.
407 specorder : list[str], optional
408 Chemical symbols in the order of LAMMPS atom types, by default None
409 force_skew : bool, optional
410 Force to write the cell as a
411 `triclinic <https://docs.lammps.org/Howto_triclinic.html>`__ box,
412 by default False
413 reduce_cell : bool, optional
414 Whether the cell shape is reduced or not, by default False
415 prismobj : Prism|None, optional
416 Prism, by default None
417 write_image_flags : bool, default False
418 If True, the image flags, i.e., in which images of the periodic
419 simulation box the atoms are, are written.
420 masses : bool, optional
421 Whether the atomic masses are written or not, by default False
422 velocities : bool, optional
423 Whether the atomic velocities are written or not, by default False
424 units : str, optional
425 `LAMMPS units <https://docs.lammps.org/units.html>`__,
426 by default 'metal'
427 bonds : bool, optional
428 Whether the bonds are written or not. Bonds can only be written
429 for atom_style='full', by default True
430 atom_style : {'atomic', 'charge', 'full'}, optional
431 `LAMMPS atom style <https://docs.lammps.org/atom_style.html>`__,
432 by default 'atomic'
433 """
435 # FIXME: We should add a check here that the encoding of the file object
436 # is actually ascii once the 'encoding' attribute of IOFormat objects
437 # starts functioning in implementation (currently it doesn't do
438 # anything).
440 if isinstance(atoms, list):
441 if len(atoms) > 1:
442 raise ValueError(
443 'Can only write one configuration to a lammps data file!'
444 )
445 atoms = atoms[0]
447 fd.write('(written by ASE)\n\n')
449 symbols = atoms.get_chemical_symbols()
450 n_atoms = len(symbols)
451 fd.write(f'{n_atoms} atoms\n')
453 if specorder is None:
454 # This way it is assured that LAMMPS atom types are always
455 # assigned predictably according to the alphabetic order
456 species = sorted(set(symbols))
457 else:
458 # To index elements in the LAMMPS data file
459 # (indices must correspond to order in the potential file)
460 species = specorder
461 n_atom_types = len(species)
462 fd.write(f'{n_atom_types} atom types\n\n')
464 bonds_in = []
465 if (bonds and (atom_style == 'full') and
466 (atoms.arrays.get('bonds') is not None)):
467 n_bonds = 0
468 n_bond_types = 1
469 for i, bondsi in enumerate(atoms.arrays['bonds']):
470 if bondsi != '_':
471 for bond in bondsi.split(','):
472 dummy1, dummy2 = bond.split('(')
473 bond_type = int(dummy2.split(')')[0])
474 at1 = int(i) + 1
475 at2 = int(dummy1) + 1
476 bonds_in.append((bond_type, at1, at2))
477 n_bonds = n_bonds + 1
478 if bond_type > n_bond_types:
479 n_bond_types = bond_type
480 fd.write(f'{n_bonds} bonds\n')
481 fd.write(f'{n_bond_types} bond types\n\n')
483 if prismobj is None:
484 prismobj = Prism(atoms.get_cell(), reduce_cell=reduce_cell)
486 # Get cell parameters and convert from ASE units to LAMMPS units
487 xhi, yhi, zhi, xy, xz, yz = convert(
488 prismobj.get_lammps_prism(), 'distance', 'ASE', units)
490 fd.write(f'0.0 {xhi:23.17g} xlo xhi\n')
491 fd.write(f'0.0 {yhi:23.17g} ylo yhi\n')
492 fd.write(f'0.0 {zhi:23.17g} zlo zhi\n')
494 if force_skew or prismobj.is_skewed():
495 fd.write(f'{xy:23.17g} {xz:23.17g} {yz:23.17g} xy xz yz\n')
496 fd.write('\n')
498 if masses:
499 _write_masses(fd, atoms, species, units)
501 # Write (unwrapped) atomic positions. If wrapping of atoms back into the
502 # cell along periodic directions is desired, this should be done manually
503 # on the Atoms object itself beforehand.
504 fd.write(f'Atoms # {atom_style}\n\n')
506 if write_image_flags:
507 scaled_positions = atoms.get_scaled_positions(wrap=False)
508 image_flags = np.floor(scaled_positions).astype(int)
510 # when `write_image_flags` is True, the positions are wrapped while the
511 # unwrapped positions can be recovered from the image flags
512 pos = prismobj.vector_to_lammps(
513 atoms.get_positions(),
514 wrap=write_image_flags,
515 )
517 if atom_style == 'atomic':
518 # Convert position from ASE units to LAMMPS units
519 pos = convert(pos, 'distance', 'ASE', units)
520 for i, r in enumerate(pos):
521 s = species.index(symbols[i]) + 1
522 line = (
523 f'{i + 1:>6} {s:>3}'
524 f' {r[0]:23.17g} {r[1]:23.17g} {r[2]:23.17g}'
525 )
526 if write_image_flags:
527 img = image_flags[i]
528 line += f' {img[0]:6d} {img[1]:6d} {img[2]:6d}'
529 line += '\n'
530 fd.write(line)
531 elif atom_style == 'charge':
532 charges = atoms.get_initial_charges()
533 # Convert position and charge from ASE units to LAMMPS units
534 pos = convert(pos, 'distance', 'ASE', units)
535 charges = convert(charges, 'charge', 'ASE', units)
536 for i, (q, r) in enumerate(zip(charges, pos)):
537 s = species.index(symbols[i]) + 1
538 line = (
539 f'{i + 1:>6} {s:>3} {q:>5}'
540 f' {r[0]:23.17g} {r[1]:23.17g} {r[2]:23.17g}'
541 )
542 if write_image_flags:
543 img = image_flags[i]
544 line += f' {img[0]:6d} {img[1]:6d} {img[2]:6d}'
545 line += '\n'
546 fd.write(line)
547 elif atom_style == 'full':
548 charges = atoms.get_initial_charges()
549 # The label 'mol-id' has apparenlty been introduced in read earlier,
550 # but so far not implemented here. Wouldn't a 'underscored' label
551 # be better, i.e. 'mol_id' or 'molecule_id'?
552 if atoms.has('mol-id'):
553 molecules = atoms.get_array('mol-id')
554 if not np.issubdtype(molecules.dtype, np.integer):
555 raise TypeError(
556 f'If "atoms" object has "mol-id" array, then '
557 f'mol-id dtype must be subtype of np.integer, and '
558 f'not {molecules.dtype!s:s}.')
559 if (len(molecules) != len(atoms)) or (molecules.ndim != 1):
560 raise TypeError(
561 'If "atoms" object has "mol-id" array, then '
562 'each atom must have exactly one mol-id.')
563 else:
564 # Assigning each atom to a distinct molecule id would seem
565 # preferableabove assigning all atoms to a single molecule
566 # id per default, as done within ase <= v 3.19.1. I.e.,
567 # molecules = np.arange(start=1, stop=len(atoms)+1,
568 # step=1, dtype=int) However, according to LAMMPS default
569 # behavior,
570 molecules = np.zeros(len(atoms), dtype=int)
571 # which is what happens if one creates new atoms within LAMMPS
572 # without explicitly taking care of the molecule id.
573 # Quote from docs at https://lammps.sandia.gov/doc/read_data.html:
574 # The molecule ID is a 2nd identifier attached to an atom.
575 # Normally, it is a number from 1 to N, identifying which
576 # molecule the atom belongs to. It can be 0 if it is a
577 # non-bonded atom or if you don't care to keep track of molecule
578 # assignments.
580 # Convert position and charge from ASE units to LAMMPS units
581 pos = convert(pos, 'distance', 'ASE', units)
582 charges = convert(charges, 'charge', 'ASE', units)
583 for i, (m, q, r) in enumerate(zip(molecules, charges, pos)):
584 s = species.index(symbols[i]) + 1
585 line = (
586 f'{i + 1:>6} {m:>3} {s:>3} {q:>5}'
587 f' {r[0]:23.17g} {r[1]:23.17g} {r[2]:23.17g}'
588 )
589 if write_image_flags:
590 img = image_flags[i]
591 line += f' {img[0]:6d} {img[1]:6d} {img[2]:6d}'
592 line += '\n'
593 fd.write(line)
594 if bonds and (atoms.arrays.get('bonds') is not None):
595 fd.write('\nBonds\n\n')
596 for i in range(n_bonds):
597 bond_type = bonds_in[i][0]
598 at1 = bonds_in[i][1]
599 at2 = bonds_in[i][2]
600 fd.write(f'{i + 1:>3} {bond_type:>3} {at1:>3} {at2:>3}\n')
601 else:
602 raise ValueError(atom_style)
604 if velocities and atoms.get_velocities() is not None:
605 fd.write('\n\nVelocities\n\n')
606 vel = prismobj.vector_to_lammps(atoms.get_velocities())
607 # Convert velocity from ASE units to LAMMPS units
608 vel = convert(vel, 'velocity', 'ASE', units)
609 for i, v in enumerate(vel):
610 fd.write(f'{i + 1:>6} {v[0]:23.17g} {v[1]:23.17g} {v[2]:23.17g}\n')
612 fd.flush()
615def _write_masses(fd, atoms: Atoms, species: list, units: str):
616 symbols_indices = atoms.symbols.indices()
617 fd.write('Masses\n\n')
618 for i, s in enumerate(species):
619 if s in symbols_indices:
620 # Find the first atom of the element `s` and extract its mass
621 # Cover by `float` to make a new object for safety
622 mass = float(atoms[symbols_indices[s][0]].mass)
623 else:
624 # Fetch from ASE data if the element `s` is not in the system
625 mass = atomic_masses[atomic_numbers[s]]
626 # Convert mass from ASE units to LAMMPS units
627 mass = convert(mass, 'mass', 'ASE', units)
628 atom_type = i + 1
629 fd.write(f'{atom_type:>6} {mass:23.17g} # {s}\n')
630 fd.write('\n')