1# lammps.py (2011/03/29) 

2# An ASE calculator for the LAMMPS classical MD code available from 

3# http://lammps.sandia.gov/ 

4# The environment variable ASE_LAMMPSRUN_COMMAND must be defined to point to the 

5# LAMMPS binary. 

6# 

7# Copyright (C) 2009 - 2011 Joerg Meyer, joerg.meyer@ch.tum.de 

8# 

9# This library is free software; you can redistribute it and/or 

10# modify it under the terms of the GNU Lesser General Public 

11# License as published by the Free Software Foundation; either 

12# version 2.1 of the License, or (at your option) any later version. 

13# 

14# This library is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 

17# Lesser General Public License for more details. 

18# 

19# You should have received a copy of the GNU Lesser General Public 

20# License along with this file; if not, write to the Free Software 

21# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 

22# USA or see <http://www.gnu.org/licenses/>. 

23 

24 

25import os 

26import shutil 

27import shlex 

28from subprocess import Popen, PIPE, TimeoutExpired 

29from threading import Thread 

30from re import compile as re_compile, IGNORECASE 

31from tempfile import mkdtemp, NamedTemporaryFile, mktemp as uns_mktemp 

32import inspect 

33import warnings 

34from typing import Dict, Any 

35import numpy as np 

36 

37from ase import Atoms 

38from ase.parallel import paropen 

39from ase.calculators.calculator import Calculator 

40from ase.calculators.calculator import all_changes 

41from ase.data import chemical_symbols 

42from ase.data import atomic_masses 

43from ase.io.lammpsdata import write_lammps_data 

44from ase.io.lammpsrun import read_lammps_dump 

45from ase.calculators.lammps import Prism 

46from ase.calculators.lammps import write_lammps_in 

47from ase.calculators.lammps import CALCULATION_END_MARK 

48from ase.calculators.lammps import convert 

49 

50__all__ = ["LAMMPS"] 

51 

52 

53class LAMMPS(Calculator): 

54 """The LAMMPS calculators object 

55 

56 files: list 

57 List of files typically containing relevant potentials for the 

58 calculation 

59 parameters: dict 

60 Dictionary of settings to be passed into the input file for calculation. 

61 specorder: list 

62 Within LAAMPS, atoms are identified by an integer value starting from 1. 

63 This variable allows the user to define the order of the indices 

64 assigned to the atoms in the calculation, with the default 

65 if not given being alphabetical 

66 keep_tmp_files: bool 

67 Retain any temporary files created. Mostly useful for debugging. 

68 tmp_dir: str 

69 path/dirname (default None -> create automatically). 

70 Explicitly control where the calculator object should create 

71 its files. Using this option implies 'keep_tmp_files' 

72 no_data_file: bool 

73 Controls whether an explicit data file will be used for feeding 

74 atom coordinates into lammps. Enable it to lessen the pressure on 

75 the (tmp) file system. THIS OPTION MIGHT BE UNRELIABLE FOR CERTAIN 

76 CORNER CASES (however, if it fails, you will notice...). 

77 keep_alive: bool 

78 When using LAMMPS as a spawned subprocess, keep the subprocess 

79 alive (but idling when unused) along with the calculator object. 

80 always_triclinic: bool 

81 Force use of a triclinic cell in LAMMPS, even if the cell is 

82 a perfect parallelepiped. 

83 

84 **Example** 

85 

86Provided that the respective potential file is in the working directory, one 

87can simply run (note that LAMMPS needs to be compiled to work with EAM 

88potentials) 

89 

90:: 

91 

92 from ase import Atom, Atoms 

93 from ase.build import bulk 

94 from ase.calculators.lammpsrun import LAMMPS 

95 

96 parameters = {'pair_style': 'eam/alloy', 

97 'pair_coeff': ['* * NiAlH_jea.eam.alloy H Ni']} 

98 

99 files = ['NiAlH_jea.eam.alloy'] 

100 

101 Ni = bulk('Ni', cubic=True) 

102 H = Atom('H', position=Ni.cell.diagonal()/2) 

103 NiH = Ni + H 

104 

105 lammps = LAMMPS(parameters=parameters, files=files) 

106 

107 NiH.calc = lammps 

108 print("Energy ", NiH.get_potential_energy()) 

109 

110(Remember you also need to set the environment variable 

111``$ASE_LAMMPSRUN_COMMAND``) 

112 

113 """ 

114 

115 name = "lammpsrun" 

116 implemented_properties = ["energy", "forces", "stress", "energies"] 

117 

118 # parameters to choose options in LAMMPSRUN 

119 ase_parameters: Dict[str, Any] = dict( 

120 specorder=None, 

121 always_triclinic=False, 

122 keep_alive=True, 

123 keep_tmp_files=False, 

124 no_data_file=False, 

125 tmp_dir=None, 

126 files=[], # usually contains potential parameters 

127 verbose=False, 

128 write_velocities=False, 

129 binary_dump=True, # bool - use binary dump files (full 

130 # precision but long long ids are casted to 

131 # double) 

132 lammps_options="-echo log -screen none -log /dev/stdout", 

133 trajectory_out=None, # file object, if is not None the 

134 # trajectory will be saved in it 

135 ) 

136 

137 # parameters forwarded to LAMMPS 

138 lammps_parameters = dict( 

139 boundary=None, # bounadry conditions styles 

140 units="metal", # str - Which units used; some potentials 

141 # require certain units 

142 atom_style="atomic", 

143 special_bonds=None, 

144 # potential informations 

145 pair_style="lj/cut 2.5", 

146 pair_coeff=["* * 1 1"], 

147 masses=None, 

148 pair_modify=None, 

149 # variables controlling the output 

150 thermo_args=[ 

151 "step", "temp", "press", "cpu", "pxx", "pyy", "pzz", 

152 "pxy", "pxz", "pyz", "ke", "pe", "etotal", "vol", "lx", 

153 "ly", "lz", "atoms", ], 

154 dump_properties=["id", "type", "x", "y", "z", "vx", "vy", 

155 "vz", "fx", "fy", "fz", ], 

156 dump_period=1, # period of system snapshot saving (in MD steps) 

157 ) 

158 

159 default_parameters = dict(ase_parameters, **lammps_parameters) 

160 

161 # legacy parameter persist, when the 'parameters' dictinary is manipulated 

162 # from the outside. All others are rested to the default value 

163 legacy_parameters = [ 

164 "specorder", 

165 "dump_period", 

166 "always_triclinic", 

167 "keep_alive", 

168 "keep_tmp_files", 

169 "tmp_dir", 

170 "parameters", 

171 "no_data_file", 

172 "files", 

173 "write_velocities", 

174 "trajectory_out", 

175 ] 

176 

177 legacy_parameters_map = {"_custom_thermo_args": "thermo_args"} 

178 

179 legacy_warn_string = "You are using an " 

180 legacy_warn_string += "old syntax to set '{}'.\n" 

181 legacy_warn_string += "Please use {}.set().".format(name.upper()) 

182 

183 def __init__(self, label="lammps", **kwargs): 

184 # "Parameters" used to be the dictionary with all parameters forwarded 

185 # to lammps. This clashes with the implementation in Calculator to 

186 # reload an old one. Trying to catch both cases to not break old 

187 # scripts. 

188 if "parameters" in kwargs: 

189 old_parameters = kwargs["parameters"] 

190 if isinstance(old_parameters, dict): 

191 warnings.warn(self.legacy_warn_string.format("parameters")) 

192 del kwargs["parameters"] 

193 else: 

194 old_parameters = None 

195 

196 Calculator.__init__(self, label=label, **kwargs) 

197 

198 if old_parameters and isinstance(old_parameters, dict): 

199 self.set(**old_parameters) 

200 

201 self.prism = None 

202 self.calls = 0 

203 self.forces = None 

204 # thermo_content contains data "written by" thermo_style. 

205 # It is a list of dictionaries, each dict (one for each line 

206 # printed by thermo_style) contains a mapping between each 

207 # custom_thermo_args-argument and the corresponding 

208 # value as printed by lammps. thermo_content will be 

209 # re-populated by the read_log method. 

210 self.thermo_content = [] 

211 

212 if self.parameters.tmp_dir is not None: 

213 # If tmp_dir is pointing somewhere, don't remove stuff! 

214 self.parameters.keep_tmp_files = True 

215 self._lmp_handle = None # To handle the lmp process 

216 

217 if self.parameters.tmp_dir is None: 

218 self.parameters.tmp_dir = mkdtemp(prefix="LAMMPS-") 

219 else: 

220 self.parameters.tmp_dir = os.path.realpath(self.parameters.tmp_dir) 

221 if not os.path.isdir(self.parameters.tmp_dir): 

222 os.mkdir(self.parameters.tmp_dir, 0o755) 

223 

224 for f in self.parameters.files: 

225 shutil.copy( 

226 f, os.path.join(self.parameters.tmp_dir, os.path.basename(f)) 

227 ) 

228 

229 def get_lammps_command(self): 

230 cmd = self.parameters.get('command') 

231 if cmd is None: 

232 envvar = 'ASE_{}_COMMAND'.format(self.name.upper()) 

233 cmd = os.environ.get(envvar) 

234 

235 if cmd is None: 

236 cmd = 'lammps' 

237 

238 opts = self.parameters.get('lammps_options') 

239 

240 if opts is not None: 

241 cmd = '{} {}'.format(cmd, opts) 

242 

243 return cmd 

244 

245 def __setattr__(self, key, value): 

246 """Catch attribute sets to emulate legacy behavior. 

247 

248 Old LAMMPSRUN allows to just override the parameters 

249 dictionary. "Modern" ase calculators can assume that default 

250 parameters are always set, overrides of the 

251 'parameters'-dictionary have to be caught and the default 

252 parameters need to be added first. A check refuses to set 

253 calculator attributes if they are unknown and set outside the 

254 '__init__' functions. 

255 """ 

256 # !TODO: remove and break somebody's code (e.g. the test examples) 

257 if ( 

258 key == "parameters" 

259 and value is not None 

260 and self.parameters is not None 

261 ): 

262 temp_dict = self.get_default_parameters() 

263 if self.parameters: 

264 for l_key in self.legacy_parameters: 

265 try: 

266 temp_dict[l_key] = self.parameters[l_key] 

267 except KeyError: 

268 pass 

269 temp_dict.update(value) 

270 value = temp_dict 

271 if key in self.legacy_parameters and key != "parameters": 

272 warnings.warn(self.legacy_warn_string.format(key)) 

273 self.set(**{key: value}) 

274 elif key in self.legacy_parameters_map: 

275 warnings.warn( 

276 self.legacy_warn_string.format( 

277 "{} for {}".format(self.legacy_parameters_map[key], key) 

278 ) 

279 ) 

280 self.set(**{self.legacy_parameters_map[key]: value}) 

281 # Catch setting none-default attributes 

282 # one test was assigning an useless Attribute, but it still worked 

283 # because the assigned object was before manipulation already handed 

284 # over to the calculator (10/2018) 

285 elif hasattr(self, key) or inspect.stack()[1][3] == "__init__": 

286 Calculator.__setattr__(self, key, value) 

287 else: 

288 raise AttributeError("Setting unknown Attribute '{}'".format(key)) 

289 

290 def __getattr__(self, key): 

291 """Corresponding getattribute-function to emulate legacy behavior. 

292 """ 

293 if key in self.legacy_parameters and key != "parameters": 

294 return self.parameters[key] 

295 if key in self.legacy_parameters_map: 

296 return self.parameters[self.legacy_parameters_map[key]] 

297 return object.__getattribute__(self, key) 

298 

299 def clean(self, force=False): 

300 

301 self._lmp_end() 

302 

303 if not self.parameters.keep_tmp_files or force: 

304 shutil.rmtree(self.parameters.tmp_dir) 

305 

306 def check_state(self, atoms, tol=1.0e-10): 

307 # Transforming the unit cell to conform to LAMMPS' convention for 

308 # orientation (c.f. https://lammps.sandia.gov/doc/Howto_triclinic.html) 

309 # results in some precision loss, so we use bit larger tolerance than 

310 # machine precision here. Note that there can also be precision loss 

311 # related to how many significant digits are specified for things in 

312 # the LAMMPS input file. 

313 return Calculator.check_state(self, atoms, tol) 

314 

315 def calculate(self, atoms=None, properties=None, system_changes=None): 

316 if properties is None: 

317 properties = self.implemented_properties 

318 if system_changes is None: 

319 system_changes = all_changes 

320 Calculator.calculate(self, atoms, properties, system_changes) 

321 self.run() 

322 

323 def _lmp_alive(self): 

324 # Return True if this calculator is currently handling a running 

325 # lammps process 

326 return self._lmp_handle and not isinstance( 

327 self._lmp_handle.poll(), int 

328 ) 

329 

330 def _lmp_end(self): 

331 # Close lammps input and wait for lammps to end. Return process 

332 # return value 

333 if self._lmp_alive(): 

334 # !TODO: handle lammps error codes 

335 try: 

336 self._lmp_handle.communicate(timeout=5) 

337 except TimeoutExpired: 

338 self._lmp_handle.kill() 

339 self._lmp_handle.communicate() 

340 err = self._lmp_handle.poll() 

341 assert err is not None 

342 return err 

343 

344 def set_missing_parameters(self): 

345 """Verify that all necessary variables are set. 

346 """ 

347 symbols = self.atoms.get_chemical_symbols() 

348 # If unspecified default to atom types in alphabetic order 

349 if not self.parameters.specorder: 

350 self.parameters.specorder = sorted(set(symbols)) 

351 

352 # !TODO: handle cases were setting masses actual lead to errors 

353 if not self.parameters.masses: 

354 self.parameters.masses = [] 

355 for type_id, specie in enumerate(self.parameters.specorder): 

356 mass = atomic_masses[chemical_symbols.index(specie)] 

357 self.parameters.masses += [ 

358 "{0:d} {1:f}".format(type_id + 1, mass) 

359 ] 

360 

361 # set boundary condtions 

362 if not self.parameters.boundary: 

363 b_str = " ".join(["fp"[int(x)] for x in self.atoms.get_pbc()]) 

364 self.parameters.boundary = b_str 

365 

366 def run(self, set_atoms=False): 

367 # !TODO: split this function 

368 """Method which explicitly runs LAMMPS.""" 

369 pbc = self.atoms.get_pbc() 

370 if all(pbc): 

371 cell = self.atoms.get_cell() 

372 elif not any(pbc): 

373 # large enough cell for non-periodic calculation - 

374 # LAMMPS shrink-wraps automatically via input command 

375 # "periodic s s s" 

376 # below 

377 cell = 2 * np.max(np.abs(self.atoms.get_positions())) * np.eye(3) 

378 else: 

379 warnings.warn( 

380 "semi-periodic ASE cell detected - translation " 

381 + "to proper LAMMPS input cell might fail" 

382 ) 

383 cell = self.atoms.get_cell() 

384 self.prism = Prism(cell) 

385 

386 self.set_missing_parameters() 

387 self.calls += 1 

388 

389 # change into subdirectory for LAMMPS calculations 

390 tempdir = self.parameters.tmp_dir 

391 

392 # setup file names for LAMMPS calculation 

393 label = "{0}{1:>06}".format(self.label, self.calls) 

394 lammps_in = uns_mktemp( 

395 prefix="in_" + label, dir=tempdir 

396 ) 

397 lammps_log = uns_mktemp( 

398 prefix="log_" + label, dir=tempdir 

399 ) 

400 lammps_trj_fd = NamedTemporaryFile( 

401 prefix="trj_" + label, 

402 suffix=(".bin" if self.parameters.binary_dump else ""), 

403 dir=tempdir, 

404 delete=(not self.parameters.keep_tmp_files), 

405 ) 

406 lammps_trj = lammps_trj_fd.name 

407 if self.parameters.no_data_file: 

408 lammps_data = None 

409 else: 

410 lammps_data_fd = NamedTemporaryFile( 

411 prefix="data_" + label, 

412 dir=tempdir, 

413 delete=(not self.parameters.keep_tmp_files), 

414 mode='w', 

415 encoding='ascii' 

416 ) 

417 write_lammps_data( 

418 lammps_data_fd, 

419 self.atoms, 

420 specorder=self.parameters.specorder, 

421 force_skew=self.parameters.always_triclinic, 

422 velocities=self.parameters.write_velocities, 

423 prismobj=self.prism, 

424 units=self.parameters.units, 

425 atom_style=self.parameters.atom_style 

426 ) 

427 lammps_data = lammps_data_fd.name 

428 lammps_data_fd.flush() 

429 

430 # see to it that LAMMPS is started 

431 if not self._lmp_alive(): 

432 command = self.get_lammps_command() 

433 # Attempt to (re)start lammps 

434 self._lmp_handle = Popen( 

435 shlex.split(command, posix=(os.name == "posix")), 

436 stdin=PIPE, 

437 stdout=PIPE, 

438 ) 

439 lmp_handle = self._lmp_handle 

440 

441 # Create thread reading lammps stdout (for reference, if requested, 

442 # also create lammps_log, although it is never used) 

443 if self.parameters.keep_tmp_files: 

444 lammps_log_fd = open(lammps_log, "wb") 

445 fd = SpecialTee(lmp_handle.stdout, lammps_log_fd) 

446 else: 

447 fd = lmp_handle.stdout 

448 thr_read_log = Thread(target=self.read_lammps_log, args=(fd,)) 

449 thr_read_log.start() 

450 

451 # write LAMMPS input (for reference, also create the file lammps_in, 

452 # although it is never used) 

453 if self.parameters.keep_tmp_files: 

454 lammps_in_fd = open(lammps_in, "wb") 

455 fd = SpecialTee(lmp_handle.stdin, lammps_in_fd) 

456 else: 

457 fd = lmp_handle.stdin 

458 write_lammps_in( 

459 lammps_in=fd, 

460 parameters=self.parameters, 

461 atoms=self.atoms, 

462 prismobj=self.prism, 

463 lammps_trj=lammps_trj, 

464 lammps_data=lammps_data, 

465 ) 

466 

467 if self.parameters.keep_tmp_files: 

468 lammps_in_fd.close() 

469 

470 # Wait for log output to be read (i.e., for LAMMPS to finish) 

471 # and close the log file if there is one 

472 thr_read_log.join() 

473 if self.parameters.keep_tmp_files: 

474 lammps_log_fd.close() 

475 

476 if not self.parameters.keep_alive: 

477 self._lmp_end() 

478 

479 exitcode = lmp_handle.poll() 

480 if exitcode and exitcode != 0: 

481 raise RuntimeError( 

482 "LAMMPS exited in {} with exit code: {}." 

483 "".format(tempdir, exitcode) 

484 ) 

485 

486 # A few sanity checks 

487 if len(self.thermo_content) == 0: 

488 raise RuntimeError("Failed to retrieve any thermo_style-output") 

489 if int(self.thermo_content[-1]["atoms"]) != len(self.atoms): 

490 # This obviously shouldn't happen, but if prism.fold_...() fails, 

491 # it could 

492 raise RuntimeError("Atoms have gone missing") 

493 

494 trj_atoms = read_lammps_dump( 

495 infileobj=lammps_trj, 

496 order=False, 

497 index=-1, 

498 prismobj=self.prism, 

499 specorder=self.parameters.specorder, 

500 ) 

501 

502 if set_atoms: 

503 self.atoms = trj_atoms.copy() 

504 

505 self.forces = trj_atoms.get_forces() 

506 # !TODO: trj_atoms is only the last snapshot of the system; Is it 

507 # desirable to save also the inbetween steps? 

508 if self.parameters.trajectory_out is not None: 

509 # !TODO: is it advisable to create here temporary atoms-objects 

510 self.trajectory_out.write(trj_atoms) 

511 

512 tc = self.thermo_content[-1] 

513 self.results["energy"] = convert( 

514 tc["pe"], "energy", self.parameters["units"], "ASE" 

515 ) 

516 self.results["free_energy"] = self.results["energy"] 

517 self.results["forces"] = self.forces.copy() 

518 stress = np.array( 

519 [-tc[i] for i in ("pxx", "pyy", "pzz", "pyz", "pxz", "pxy")] 

520 ) 

521 

522 # We need to apply the Lammps rotation stuff to the stress: 

523 xx, yy, zz, yz, xz, xy = stress 

524 stress_tensor = np.array([[xx, xy, xz], 

525 [xy, yy, yz], 

526 [xz, yz, zz]]) 

527 R = self.prism.rot_mat 

528 stress_atoms = np.dot(R, stress_tensor) 

529 stress_atoms = np.dot(stress_atoms, R.T) 

530 stress_atoms = stress_atoms[[0, 1, 2, 1, 0, 0], 

531 [0, 1, 2, 2, 2, 1]] 

532 stress = stress_atoms 

533 

534 self.results["stress"] = convert( 

535 stress, "pressure", self.parameters["units"], "ASE" 

536 ) 

537 

538 lammps_trj_fd.close() 

539 if not self.parameters.no_data_file: 

540 lammps_data_fd.close() 

541 

542 def __enter__(self): 

543 return self 

544 

545 def __exit__(self, *args): 

546 self._lmp_end() 

547 

548 def read_lammps_log(self, lammps_log=None): 

549 # !TODO: somehow communicate 'thermo_content' explicitly 

550 """Method which reads a LAMMPS output log file.""" 

551 

552 if lammps_log is None: 

553 lammps_log = self.label + ".log" 

554 

555 if isinstance(lammps_log, str): 

556 fileobj = paropen(lammps_log, "wb") 

557 close_log_file = True 

558 else: 

559 # Expect lammps_in to be a file-like object 

560 fileobj = lammps_log 

561 close_log_file = False 

562 

563 # read_log depends on that the first (three) thermo_style custom args 

564 # can be capitalized and matched against the log output. I.e. 

565 # don't use e.g. 'ke' or 'cpu' which are labeled KinEng and CPU. 

566 _custom_thermo_mark = " ".join( 

567 [x.capitalize() for x in self.parameters.thermo_args[0:3]] 

568 ) 

569 

570 # !TODO: regex-magic necessary? 

571 # Match something which can be converted to a float 

572 f_re = r"([+-]?(?:(?:\d+(?:\.\d*)?|\.\d+)(?:e[+-]?\d+)?|nan|inf))" 

573 n_args = len(self.parameters["thermo_args"]) 

574 # Create a re matching exactly N white space separated floatish things 

575 _custom_thermo_re = re_compile( 

576 r"^\s*" + r"\s+".join([f_re] * n_args) + r"\s*$", flags=IGNORECASE 

577 ) 

578 

579 thermo_content = [] 

580 line = fileobj.readline().decode("utf-8") 

581 while line and line.strip() != CALCULATION_END_MARK: 

582 # check error 

583 if 'ERROR:' in line: 

584 if close_log_file: 

585 fileobj.close() 

586 raise RuntimeError(f'LAMMPS exits with error message: {line}') 

587 

588 # get thermo output 

589 if line.startswith(_custom_thermo_mark): 

590 bool_match = True 

591 while bool_match: 

592 line = fileobj.readline().decode("utf-8") 

593 bool_match = _custom_thermo_re.match(line) 

594 if bool_match: 

595 # create a dictionary between each of the 

596 # thermo_style args and it's corresponding value 

597 thermo_content.append( 

598 dict( 

599 zip( 

600 self.parameters.thermo_args, 

601 map(float, bool_match.groups()), 

602 ) 

603 ) 

604 ) 

605 else: 

606 line = fileobj.readline().decode("utf-8") 

607 

608 if close_log_file: 

609 fileobj.close() 

610 

611 self.thermo_content = thermo_content 

612 

613 

614class SpecialTee: 

615 """A special purpose, with limited applicability, tee-like thing. 

616 

617 A subset of stuff read from, or written to, orig_fd, 

618 is also written to out_fd. 

619 It is used by the lammps calculator for creating file-logs of stuff 

620 read from, or written to, stdin and stdout, respectively. 

621 """ 

622 

623 def __init__(self, orig_fd, out_fd): 

624 self._orig_fd = orig_fd 

625 self._out_fd = out_fd 

626 self.name = orig_fd.name 

627 

628 def write(self, data): 

629 self._orig_fd.write(data) 

630 self._out_fd.write(data) 

631 self.flush() 

632 

633 def read(self, *args, **kwargs): 

634 data = self._orig_fd.read(*args, **kwargs) 

635 self._out_fd.write(data) 

636 return data 

637 

638 def readline(self, *args, **kwargs): 

639 data = self._orig_fd.readline(*args, **kwargs) 

640 self._out_fd.write(data) 

641 return data 

642 

643 def readlines(self, *args, **kwargs): 

644 data = self._orig_fd.readlines(*args, **kwargs) 

645 self._out_fd.write("".join(data)) 

646 return data 

647 

648 def flush(self): 

649 self._orig_fd.flush() 

650 self._out_fd.flush() 

651 

652 

653if __name__ == "__main__": 

654 pair_style = "eam" 

655 Pd_eam_file = "Pd_u3.eam" 

656 pair_coeff = ["* * " + Pd_eam_file] 

657 parameters = {"pair_style": pair_style, "pair_coeff": pair_coeff} 

658 my_files = [Pd_eam_file] 

659 calc = LAMMPS(parameters=parameters, files=my_files) 

660 a0 = 3.93 

661 b0 = a0 / 2.0 

662 

663 bulk = Atoms( 

664 ["Pd"] * 4, 

665 positions=[(0, 0, 0), (b0, b0, 0), (b0, 0, b0), (0, b0, b0)], 

666 cell=[a0] * 3, 

667 pbc=True, 

668 ) 

669 # test get_forces 

670 print("forces for a = {0}".format(a0)) 

671 print(calc.get_forces(bulk)) 

672 # single points for various lattice constants 

673 bulk.calc = calc 

674 for i in range(-5, 5, 1): 

675 a = a0 * (1 + i / 100.0) 

676 bulk.set_cell([a] * 3) 

677 print("a : {0} , total energy : {1}".format( 

678 a, bulk.get_potential_energy())) 

679 

680 calc.clean()