 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1'''Constant pressure/stress and temperature dynamics.

3Combined Nose-Hoover and Parrinello-Rahman dynamics, creating an NPT

4(or N,stress,T) ensemble.

6The method is the one proposed by Melchionna et al.  and later

7modified by Melchionna . The differential equations are integrated

8using a centered difference method .

10 1. S. Melchionna, G. Ciccotti and B. L. Holian, "Hoover NPT dynamics

11 for systems varying in shape and size", Molecular Physics 78, p. 533

12 (1993).

14 2. S. Melchionna, "Constrained systems and statistical distribution",

15 Physical Review E, 61, p. 6165 (2000).

17 3. B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover,

18 "Time-reversible equilibrium and nonequilibrium isothermal-isobaric

19 simulations with centered-difference Stoermer algorithms.", Physical

20 Review A, 41, p. 4552 (1990).

21'''

23import sys

24import weakref

26import numpy as np

28from ase.md.md import MolecularDynamics

29from ase import units

31linalg = np.linalg

33# Delayed imports: If the trajectory object is reading a special ASAP version

34# of HooverNPT, that class is imported from Asap.Dynamics.NPTDynamics.

37class NPT(MolecularDynamics):

39 classname = "NPT" # Used by the trajectory.

40 _npt_version = 2 # Version number, used for Asap compatibility.

42 def __init__(self, atoms,

43 timestep, temperature=None, externalstress=None,

44 ttime=None, pfactor=None,

45 *, temperature_K=None,

47 append_trajectory=False):

48 '''Constant pressure/stress and temperature dynamics.

50 Combined Nose-Hoover and Parrinello-Rahman dynamics, creating an

51 NPT (or N,stress,T) ensemble.

53 The method is the one proposed by Melchionna et al.  and later

54 modified by Melchionna . The differential equations are integrated

57 The dynamics object is called with the following parameters:

59 atoms: Atoms object

60 The list of atoms.

62 timestep: float

63 The timestep in units matching eV, Å, u.

65 temperature: float (deprecated)

66 The desired temperature in eV.

68 temperature_K: float

69 The desired temperature in K.

71 externalstress: float or nparray

72 The external stress in eV/A^3. Either a symmetric

73 3x3 tensor, a 6-vector representing the same, or a

74 scalar representing the pressure. Note that the

75 stress is positive in tension whereas the pressure is

76 positive in compression: giving a scalar p is

77 equivalent to giving the tensor (-p, -p, -p, 0, 0, 0).

79 ttime: float

80 Characteristic timescale of the thermostat, in ASE internal units

81 Set to None to disable the thermostat.

83 WARNING: Not specifying ttime sets it to None, disabling the

84 thermostat.

86 pfactor: float

87 A constant in the barostat differential equation. If

88 a characteristic barostat timescale of ptime is

89 desired, set pfactor to ptime^2 * B

90 (where ptime is in units matching

91 eV, Å, u; and B is the Bulk Modulus, given in eV/Å^3).

92 Set to None to disable the barostat.

93 Typical metallic bulk moduli are of the order of

94 100 GPa or 0.6 eV/A^3.

96 WARNING: Not specifying pfactor sets it to None, disabling the

97 barostat.

99 mask: None or 3-tuple or 3x3 nparray (optional)

100 Optional argument. A tuple of three integers (0 or 1),

101 indicating if the system can change size along the

102 three Cartesian axes. Set to (1,1,1) or None to allow

103 a fully flexible computational box. Set to (1,1,0)

104 to disallow elongations along the z-axis etc.

105 mask may also be specified as a symmetric 3x3 array

106 indicating which strain values may change.

108 Useful parameter values:

110 * The same timestep can be used as in Verlet dynamics, i.e. 5 fs is fine

111 for bulk copper.

113 * The ttime and pfactor are quite critical, too small values may

114 cause instabilites and/or wrong fluctuations in T / p. Too

115 large values cause an oscillation which is slow to die. Good

116 values for the characteristic times seem to be 25 fs for ttime,

117 and 75 fs for ptime (used to calculate pfactor), at least for

118 bulk copper with 15000-200000 atoms. But this is not well

119 tested, it is IMPORTANT to monitor the temperature and

120 stress/pressure fluctuations.

123 References:

125 1) S. Melchionna, G. Ciccotti and B. L. Holian, Molecular

126 Physics 78, p. 533 (1993).

128 2) S. Melchionna, Physical

129 Review E 61, p. 6165 (2000).

131 3) B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover,

132 Physical Review A 41, p. 4552 (1990).

134 4) F. D. Di Tolla and M. Ronchetti, Physical

135 Review E 48, p. 1726 (1993).

137 '''

139 MolecularDynamics.__init__(self, atoms, timestep, trajectory,

141 append_trajectory=append_trajectory)

142 # self.atoms = atoms

143 # self.timestep = timestep

144 if externalstress is None and pfactor is not None:

145 raise TypeError("Missing 'externalstress' argument.")

146 self.zero_center_of_mass_momentum(verbose=1)

147 self.temperature = units.kB * self._process_temperature(

148 temperature, temperature_K, 'eV')

149 self.set_stress(externalstress)

151 self.eta = np.zeros((3, 3), float)

152 self.zeta = 0.0

153 self.zeta_integrated = 0.0

154 self.initialized = 0

155 self.ttime = ttime

156 self.pfactor_given = pfactor

157 self._calculateconstants()

158 self.timeelapsed = 0.0

159 self.frac_traceless = 1

161 def set_temperature(self, temperature=None, *, temperature_K=None):

162 """Set the temperature.

164 Parameters:

166 temperature: float (deprecated)

167 The new temperature in eV. Deprecated, use ``temperature_K``.

169 temperature_K: float (keyword-only argument)

170 The new temperature, in K.

171 """

172 self.temperature = units.kB * self._process_temperature(

173 temperature, temperature_K, 'eV')

174 self._calculateconstants()

176 def set_stress(self, stress):

177 """Set the applied stress.

179 Must be a symmetric 3x3 tensor, a 6-vector representing a symmetric

180 3x3 tensor, or a number representing the pressure.

182 Use with care, it is better to set the correct stress when creating

183 the object.

184 """

186 if np.isscalar(stress):

187 stress = np.array([-stress, -stress, -stress, 0.0, 0.0, 0.0])

188 else:

189 stress = np.array(stress)

190 if stress.shape == (3, 3):

191 if not self._issymmetric(stress):

192 raise ValueError(

193 "The external stress must be a symmetric tensor.")

194 stress = np.array((stress[0, 0], stress[1, 1],

195 stress[2, 2], stress[1, 2],

196 stress[0, 2], stress[0, 1]))

197 elif stress.shape != (6,):

198 raise ValueError("The external stress has the wrong shape.")

199 self.externalstress = stress

202 """Set the mask indicating dynamic elements of the computational box.

204 If set to None, all elements may change. If set to a 3-vector

205 of ones and zeros, elements which are zero specify directions

206 along which the size of the computational box cannot change.

207 For example, if mask = (1,1,0) the length of the system along

208 the z-axis cannot change, although xz and yz shear is still

209 possible. May also be specified as a symmetric 3x3 array indicating

210 which strain values may change.

212 Use with care, as you may "freeze in" a fluctuation in the strain rate.

213 """

219 raise RuntimeError('The mask has the wrong shape ' +

220 '(must be a 3-vector or 3x3 matrix)')

221 else:

226 else:

229 def set_fraction_traceless(self, fracTraceless):

230 """set what fraction of the traceless part of the force

231 on eta is kept.

233 By setting this to zero, the volume may change but the shape may not.

234 """

235 self.frac_traceless = fracTraceless

237 def get_strain_rate(self):

238 """Get the strain rate as an upper-triangular 3x3 matrix.

240 This includes the fluctuations in the shape of the computational box.

242 """

243 return np.array(self.eta, copy=1)

245 def set_strain_rate(self, rate):

246 """Set the strain rate. Must be an upper triangular 3x3 matrix.

248 If you set a strain rate along a direction that is "masked out"

249 (see ``set_mask``), the strain rate along that direction will be

250 maintained constantly.

251 """

252 if not (rate.shape == (3, 3) and self._isuppertriangular(rate)):

253 raise ValueError("Strain rate must be an upper triangular matrix.")

254 self.eta = rate

255 if self.initialized:

256 # Recalculate h_past and eta_past so they match the current value.

257 self._initialize_eta_h()

259 def get_time(self):

260 "Get the elapsed time."

261 return self.timeelapsed

263 def run(self, steps):

264 """Perform a number of time steps."""

265 if not self.initialized:

266 self.initialize()

267 else:

268 if self.have_the_atoms_been_changed():

269 raise NotImplementedError(

270 "You have modified the atoms since the last timestep.")

272 for i in range(steps):

273 self.step()

274 self.nsteps += 1

275 self.call_observers()

277 def have_the_atoms_been_changed(self):

278 "Checks if the user has modified the positions or momenta of the atoms"

279 limit = 1e-10

280 h = self._getbox()

281 if max(abs((h - self.h).ravel())) > limit:

282 self._warning("The computational box has been modified.")

283 return 1

284 expected_r = np.dot(self.q + 0.5, h)

285 err = max(abs((expected_r - self.atoms.get_positions()).ravel()))

286 if err > limit:

287 self._warning("The atomic positions have been modified: " +

288 str(err))

289 return 1

290 return 0

292 def step(self):

293 """Perform a single time step.

295 Assumes that the forces and stresses are up to date, and that

296 the positions and momenta have not been changed since last

297 timestep.

298 """

300 # Assumes the following variables are OK

301 # q_past, q, q_future, p, eta, eta_past, zeta, zeta_past, h, h_past

302 #

303 # q corresponds to the current positions

304 # p must be equal to self.atoms.GetCartesianMomenta()

305 # h must be equal to self.atoms.GetUnitCell()

306 #

307 # print "Making a timestep"

308 dt = self.dt

309 h_future = self.h_past + 2 * dt * np.dot(self.h, self.eta)

310 if self.pfactor_given is None:

311 deltaeta = np.zeros(6, float)

312 else:

313 stress = self.stresscalculator()

314 deltaeta = -2 * dt * (self.pfact * linalg.det(self.h) *

315 (stress - self.externalstress))

317 if self.frac_traceless == 1:

318 eta_future = self.eta_past + self.mask * \

319 self._makeuppertriangular(deltaeta)

320 else:

321 trace_part, traceless_part = self._separatetrace(

322 self._makeuppertriangular(deltaeta))

323 eta_future = (self.eta_past + trace_part +

324 self.frac_traceless * traceless_part)

326 deltazeta = 2 * dt * self.tfact * (self.atoms.get_kinetic_energy() -

327 self.desiredEkin)

328 zeta_future = self.zeta_past + deltazeta

330 self.timeelapsed += dt

331 self.h_past = self.h

332 self.h = h_future

333 self.inv_h = linalg.inv(self.h)

334 self.q_past = self.q

335 self.q = self.q_future

336 self._setbox_and_positions(self.h, self.q)

337 self.eta_past = self.eta

338 self.eta = eta_future

339 self.zeta_past = self.zeta

340 self.zeta = zeta_future

341 self._synchronize() # for parallel simulations.

342 self.zeta_integrated += dt * self.zeta

343 force = self.forcecalculator()

344 self._calculate_q_future(force)

345 self.atoms.set_momenta(

346 np.dot(self.q_future - self.q_past, self.h / (2 * dt)) *

347 self._getmasses())

348 # self.stresscalculator()

350 def forcecalculator(self):

351 return self.atoms.get_forces(md=True)

353 def stresscalculator(self):

354 return self.atoms.get_stress(include_ideal_gas=True)

356 def initialize(self):

357 """Initialize the dynamics.

359 The dynamics requires positions etc for the two last times to

360 do a timestep, so the algorithm is not self-starting. This

361 method performs a 'backwards' timestep to generate a

362 configuration before the current.

364 This is called automatically the first time ``run()`` is called.

365 """

366 # print "Initializing the NPT dynamics."

367 dt = self.dt

368 atoms = self.atoms

369 self.h = self._getbox()

370 if not self._isuppertriangular(self.h):

371 print("I am", self)

372 print("self.h:")

373 print(self.h)

374 print("Min:", min((self.h[1, 0], self.h[2, 0], self.h[2, 1])))

375 print("Max:", max((self.h[1, 0], self.h[2, 0], self.h[2, 1])))

376 raise NotImplementedError(

377 "Can (so far) only operate on lists of atoms where the "

378 "computational box is an upper triangular matrix.")

379 self.inv_h = linalg.inv(self.h)

380 # The contents of the q arrays should migrate in parallel simulations.

381 # self._make_special_q_arrays()

382 self.q = np.dot(self.atoms.get_positions(), self.inv_h) - 0.5

383 # zeta and eta were set in __init__

384 self._initialize_eta_h()

385 deltazeta = dt * self.tfact * (atoms.get_kinetic_energy() -

386 self.desiredEkin)

387 self.zeta_past = self.zeta - deltazeta

388 self._calculate_q_past_and_future()

389 self.initialized = 1

391 def get_gibbs_free_energy(self):

392 """Return the Gibb's free energy, which is supposed to be conserved.

394 Requires that the energies of the atoms are up to date.

396 This is mainly intended as a diagnostic tool. If called before the

397 first timestep, Initialize will be called.

398 """

399 if not self.initialized:

400 self.initialize()

401 n = self._getnatoms()

402 # tretaTeta = sum(diagonal(matrixmultiply(transpose(self.eta),

403 # self.eta)))

404 contractedeta = np.sum((self.eta * self.eta).ravel())

405 gibbs = (self.atoms.get_potential_energy() +

406 self.atoms.get_kinetic_energy()

407 - np.sum(self.externalstress[0:3]) * linalg.det(self.h) / 3.0)

408 if self.ttime is not None:

409 gibbs += (1.5 * n * self.temperature *

410 (self.ttime * self.zeta)**2 +

411 3 * self.temperature * (n - 1) * self.zeta_integrated)

412 else:

413 assert self.zeta == 0.0

414 if self.pfactor_given is not None:

415 gibbs += 0.5 / self.pfact * contractedeta

416 else:

417 assert contractedeta == 0.0

418 return gibbs

420 def get_center_of_mass_momentum(self):

421 "Get the center of mass momentum."

422 return self.atoms.get_momenta().sum(0)

424 def zero_center_of_mass_momentum(self, verbose=0):

425 "Set the center of mass momentum to zero."

426 cm = self.get_center_of_mass_momentum()

427 abscm = np.sqrt(np.sum(cm * cm))

428 if verbose and abscm > 1e-4:

429 self._warning(

430 self.classname +

431 ": Setting the center-of-mass momentum to zero "

432 "(was %.6g %.6g %.6g)" % tuple(cm))

433 self.atoms.set_momenta(self.atoms.get_momenta() -

434 cm / self._getnatoms())

436 def attach_atoms(self, atoms):

437 """Assign atoms to a restored dynamics object.

439 This function must be called to set the atoms immediately after the

440 dynamics object has been read from a trajectory.

441 """

442 try:

443 self.atoms

444 except AttributeError:

445 pass

446 else:

447 raise RuntimeError("Cannot call attach_atoms on a dynamics "

449 MolecularDynamics.__init__(self, atoms, self.dt)

450 limit = 1e-6

451 h = self._getbox()

452 if max(abs((h - self.h).ravel())) > limit:

453 raise RuntimeError(

454 "The unit cell of the atoms does not match "

455 "the unit cell stored in the file.")

456 self.inv_h = linalg.inv(self.h)

457 self.q = np.dot(self.atoms.get_positions(), self.inv_h) - 0.5

458 self._calculate_q_past_and_future()

459 self.initialized = 1

461 def attach(self, function, interval=1, *args, **kwargs):

462 """Attach callback function or trajectory.

464 At every *interval* steps, call *function* with arguments

465 *args* and keyword arguments *kwargs*.

467 If *function* is a trajectory object, its write() method is

468 attached, but if *function* is a BundleTrajectory (or another

469 trajectory supporting set_extra_data(), said method is first

470 used to instruct the trajectory to also save internal

471 data from the NPT dynamics object.

472 """

473 if hasattr(function, "set_extra_data"):

474 # We are attaching a BundleTrajectory or similar

475 function.set_extra_data("npt_init",

476 WeakMethodWrapper(self, "get_init_data"),

477 once=True)

478 function.set_extra_data("npt_dynamics",

479 WeakMethodWrapper(self, "get_data"))

480 MolecularDynamics.attach(self, function, interval, *args, **kwargs)

482 def get_init_data(self):

483 "Return the data needed to initialize a new NPT dynamics."

484 return {'dt': self.dt,

485 'temperature': self.temperature,

486 'desiredEkin': self.desiredEkin,

487 'externalstress': self.externalstress,

489 'ttime': self.ttime,

490 'tfact': self.tfact,

491 'pfactor_given': self.pfactor_given,

492 'pfact': self.pfact,

493 'frac_traceless': self.frac_traceless}

495 def get_data(self):

496 "Return data needed to restore the state."

497 return {'eta': self.eta,

498 'eta_past': self.eta_past,

499 'zeta': self.zeta,

500 'zeta_past': self.zeta_past,

501 'zeta_integrated': self.zeta_integrated,

502 'h': self.h,

503 'h_past': self.h_past,

504 'timeelapsed': self.timeelapsed}

506 @classmethod

507 def read_from_trajectory(cls, trajectory, frame=-1, atoms=None):

508 """Read dynamics and atoms from trajectory (Class method).

510 Simultaneously reads the atoms and the dynamics from a BundleTrajectory,

511 including the internal data of the NPT dynamics object (automatically

512 saved when attaching a BundleTrajectory to an NPT object).

514 Arguments:

516 trajectory

517 The filename or an open BundleTrajectory object.

519 frame (optional)

520 Which frame to read. Default: the last.

522 atoms (optional, internal use only)

523 Pre-read atoms. Do not use.

524 """

525 if isinstance(trajectory, str):

526 if trajectory.endswith('/'):

527 trajectory = trajectory[:-1]

528 if trajectory.endswith('.bundle'):

529 from ase.io.bundletrajectory import BundleTrajectory

530 trajectory = BundleTrajectory(trajectory)

531 else:

532 raise ValueError(

533 f"Cannot open '{trajectory}': unsupported file format")

534 # trajectory is now a BundleTrajectory object (or compatible)

535 if atoms is None:

536 atoms = trajectory[frame]

539 dyn = cls(atoms, timestep=init_data['dt'],

540 temperature=init_data['temperature'],

541 externalstress=init_data['externalstress'],

542 ttime=init_data['ttime'],

543 pfactor=init_data['pfactor_given'],

545 dyn.desiredEkin = init_data['desiredEkin']

546 dyn.tfact = init_data['tfact']

547 dyn.pfact = init_data['pfact']

548 dyn.frac_traceless = init_data['frac_traceless']

549 for k, v in frame_data.items():

550 setattr(dyn, k, v)

551 return (dyn, atoms)

553 def _getbox(self):

554 "Get the computational box."

555 return self.atoms.get_cell()

557 def _getmasses(self):

558 "Get the masses as an Nx1 array."

559 return np.reshape(self.atoms.get_masses(), (-1, 1))

561 def _separatetrace(self, mat):

562 """return two matrices, one proportional to the identity

563 the other traceless, which sum to the given matrix

564 """

565 tracePart = ((mat + mat + mat) / 3.) * np.identity(3)

566 return tracePart, mat - tracePart

568 # A number of convenient helper methods

569 def _warning(self, text):

570 "Emit a warning."

571 sys.stderr.write("WARNING: " + text + "\n")

572 sys.stderr.flush()

574 def _calculate_q_future(self, force):

575 "Calculate future q. Needed in Timestep and Initialization."

576 dt = self.dt

577 id3 = np.identity(3)

578 alpha = (dt * dt) * np.dot(force / self._getmasses(),

579 self.inv_h)

580 beta = dt * np.dot(self.h, np.dot(self.eta + 0.5 * self.zeta * id3,

581 self.inv_h))

582 inv_b = linalg.inv(beta + id3)

583 self.q_future = np.dot(2 * self.q +

584 np.dot(self.q_past, beta - id3) + alpha,

585 inv_b)

587 def _calculate_q_past_and_future(self):

588 def ekin(p, m=self.atoms.get_masses()):

589 p2 = np.sum(p * p, -1)

590 return 0.5 * np.sum(p2 / m) / len(m)

591 p0 = self.atoms.get_momenta()

592 m = self._getmasses()

593 p = np.array(p0, copy=1)

594 dt = self.dt

595 for i in range(2):

596 self.q_past = self.q - dt * np.dot(p / m, self.inv_h)

597 self._calculate_q_future(self.atoms.get_forces(md=True))

598 p = np.dot(self.q_future - self.q_past, self.h / (2 * dt)) * m

599 e = ekin(p)

600 if e < 1e-5:

601 # The kinetic energy and momenta are virtually zero

602 return

603 p = (p0 - p) + p0

605 def _initialize_eta_h(self):

606 self.h_past = self.h - self.dt * np.dot(self.h, self.eta)

607 if self.pfactor_given is None:

608 deltaeta = np.zeros(6, float)

609 else:

610 deltaeta = (-self.dt * self.pfact * linalg.det(self.h)

611 * (self.stresscalculator() - self.externalstress))

612 if self.frac_traceless == 1:

613 self.eta_past = self.eta - self.mask * \

614 self._makeuppertriangular(deltaeta)

615 else:

616 trace_part, traceless_part = self._separatetrace(

617 self._makeuppertriangular(deltaeta))

618 self.eta_past = (self.eta - trace_part -

619 self.frac_traceless * traceless_part)

621 def _makeuppertriangular(self, sixvector):

622 "Make an upper triangular matrix from a 6-vector."

623 return np.array(((sixvector, sixvector, sixvector),

624 (0, sixvector, sixvector),

625 (0, 0, sixvector)))

627 @staticmethod

628 def _isuppertriangular(m) -> bool:

629 "Check that a matrix is on upper triangular form."

630 return m[1, 0] == m[2, 0] == m[2, 1] == 0.0

632 def _calculateconstants(self):

633 """(Re)calculate some constants when pfactor,

634 ttime or temperature have been changed."""

636 n = self._getnatoms()

637 if self.ttime is None:

638 self.tfact = 0.0

639 else:

640 self.tfact = 2.0 / (3 * n * self.temperature *

641 self.ttime * self.ttime)

642 if self.pfactor_given is None:

643 self.pfact = 0.0

644 else:

645 self.pfact = 1.0 / (self.pfactor_given * linalg.det(self._getbox()))

646 # self.pfact = 1.0/(n * self.temperature * self.ptime * self.ptime)

647 self.desiredEkin = 1.5 * (n - 1) * self.temperature

649 def _setbox_and_positions(self, h, q):

650 """Set the computational box and the positions."""

651 self.atoms.set_cell(h)

652 r = np.dot(q + 0.5, h)

653 self.atoms.set_positions(r)

655 # A few helper methods, which have been placed in separate methods

656 # so they can be replaced in the parallel version.

657 def _synchronize(self):

658 """Synchronize eta, h and zeta on all processors.

660 In a parallel simulation, eta, h and zeta are communicated

661 from the master to all slaves, to prevent numerical noise from

662 causing them to diverge.

664 In a serial simulation, do nothing.

665 """

666 pass # This is a serial simulation object. Do nothing.

668 def _getnatoms(self):

669 """Get the number of atoms.

671 In a parallel simulation, this is the total number of atoms on all

672 processors.

673 """

674 return len(self.atoms)

676 def _make_special_q_arrays(self):

677 """Make the arrays used to store data about the atoms.

679 In a parallel simulation, these are migrating arrays. In a

680 serial simulation they are ordinary Numeric arrays.

681 """

682 natoms = len(self.atoms)

683 self.q = np.zeros((natoms, 3), float)

684 self.q_past = np.zeros((natoms, 3), float)

685 self.q_future = np.zeros((natoms, 3), float)

688class WeakMethodWrapper:

689 """A weak reference to a method.

691 Create an object storing a weak reference to an instance and

692 the name of the method to call. When called, calls the method.

694 Just storing a weak reference to a bound method would not work,

695 as the bound method object would go away immediately.

696 """

698 def __init__(self, obj, method):

699 self.obj = weakref.proxy(obj)

700 self.method = method

702 def __call__(self, *args, **kwargs):

703 m = getattr(self.obj, self.method)

704 return m(*args, **kwargs)