1"""Logging for molecular dynamics.""" 

2 

3import weakref 

4 

5import ase.units as units 

6from ase.parallel import world 

7from ase.utils import IOContext 

8 

9 

10class MDLogger(IOContext): 

11 """Class for logging molecular dynamics simulations. 

12 

13 Parameters: 

14 dyn: The dynamics. Only a weak reference is kept. 

15 

16 atoms: The atoms. 

17 

18 logfile: File name or open file, "-" meaning standard output. 

19 

20 stress=False: Include stress in log. 

21 

22 peratom=False: Write energies per atom. 

23 

24 mode="a": How the file is opened if logfile is a filename. 

25 """ 

26 

27 def __init__(self, dyn, atoms, logfile, header=True, stress=False, 

28 peratom=False, mode="a"): 

29 if hasattr(dyn, "get_time"): 

30 self.dyn = weakref.proxy(dyn) 

31 else: 

32 self.dyn = None 

33 self.atoms = atoms 

34 global_natoms = atoms.get_global_number_of_atoms() 

35 self.logfile = self.openfile(logfile, comm=world, mode='a') 

36 self.stress = stress 

37 self.peratom = peratom 

38 if self.dyn is not None: 

39 self.hdr = "%-9s " % ("Time[ps]",) 

40 self.fmt = "%-10.4f " 

41 else: 

42 self.hdr = "" 

43 self.fmt = "" 

44 if self.peratom: 

45 self.hdr += "%12s %12s %12s %6s" % ("Etot/N[eV]", "Epot/N[eV]", 

46 "Ekin/N[eV]", "T[K]") 

47 self.fmt += "%12.4f %12.4f %12.4f %6.1f" 

48 else: 

49 self.hdr += "%12s %12s %12s %6s" % ("Etot[eV]", "Epot[eV]", 

50 "Ekin[eV]", "T[K]") 

51 # Choose a sensible number of decimals 

52 if global_natoms <= 100: 

53 digits = 4 

54 elif global_natoms <= 1000: 

55 digits = 3 

56 elif global_natoms <= 10000: 

57 digits = 2 

58 else: 

59 digits = 1 

60 self.fmt += 3 * ("%%12.%df " % (digits,)) + " %6.1f" 

61 if self.stress: 

62 self.hdr += (' ---------------------- stress [GPa] ' 

63 '-----------------------') 

64 self.fmt += 6 * " %10.3f" 

65 self.fmt += "\n" 

66 if header: 

67 self.logfile.write(self.hdr + "\n") 

68 

69 def __del__(self): 

70 self.close() 

71 

72 def __call__(self): 

73 epot = self.atoms.get_potential_energy() 

74 ekin = self.atoms.get_kinetic_energy() 

75 temp = self.atoms.get_temperature() 

76 global_natoms = self.atoms.get_global_number_of_atoms() 

77 if self.peratom: 

78 epot /= global_natoms 

79 ekin /= global_natoms 

80 if self.dyn is not None: 

81 t = self.dyn.get_time() / (1000 * units.fs) 

82 dat = (t,) 

83 else: 

84 dat = () 

85 dat += (epot + ekin, epot, ekin, temp) 

86 if self.stress: 

87 dat += tuple(self.atoms.get_stress( 

88 include_ideal_gas=True) / units.GPa) 

89 self.logfile.write(self.fmt % dat) 

90 self.logfile.flush()