|
| 1 | +# SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. |
| 2 | +# SPDX-License-Identifier: BSD-3-Clause |
| 3 | +# |
| 4 | +# Redistribution and use in source and binary forms, with or without |
| 5 | +# modification, are permitted provided that the following conditions are met: |
| 6 | +# |
| 7 | +# 1. Redistributions of source code must retain the above copyright notice, |
| 8 | +# this list of conditions and the following disclaimer. |
| 9 | +# |
| 10 | +# 2. Redistributions in binary form must reproduce the above copyright notice, |
| 11 | +# this list of conditions and the following disclaimer in the documentation |
| 12 | +# and/or other materials provided with the distribution. |
| 13 | +# |
| 14 | +# 3. Neither the name of the copyright holder nor the names of its contributors |
| 15 | +# may be used to endorse or promote products derived from this software |
| 16 | +# without specific prior written permission. |
| 17 | +# |
| 18 | +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
| 19 | +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| 20 | +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| 21 | +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE |
| 22 | +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
| 23 | +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
| 24 | +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
| 25 | +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
| 26 | +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
| 27 | +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| 28 | +# POSSIBILITY OF SUCH DAMAGE. |
| 29 | +""" |
| 30 | +NVT Molecular Dynamics with TensorNet and nvalchemi-toolkit |
| 31 | +============================================================ |
| 32 | +
|
| 33 | +This example runs NVT Langevin molecular dynamics on a periodic LiFePO4 |
| 34 | +crystal using a pretrained MatGL TensorNet potential wrapped for the |
| 35 | +nvalchemi-toolkit dynamics engine. |
| 36 | +
|
| 37 | +The workflow: |
| 38 | +
|
| 39 | +1. Load a pretrained TensorNet potential from MatGL. |
| 40 | +2. Wrap it with :class:`~matgl.ext.alchmtk.TensorNetWrapper`. |
| 41 | +3. Build a periodic LiFePO4 structure using pymatgen. |
| 42 | +4. Convert to :class:`~nvalchemi.data.AtomicData` via ``from_structure``. |
| 43 | +5. Run 1000 NVT Langevin steps at 300 K. |
| 44 | +6. Compute final temperature from kinetic energy. |
| 45 | +
|
| 46 | +Requirements:: |
| 47 | +
|
| 48 | + pip install matgl[alchmtk] |
| 49 | +
|
| 50 | +""" |
| 51 | + |
| 52 | +from __future__ import annotations |
| 53 | + |
| 54 | +import torch |
| 55 | +from nvalchemi.data import AtomicData, Batch |
| 56 | +from nvalchemi.dynamics import NVTLangevin |
| 57 | +from nvalchemi.dynamics.base import DynamicsStage |
| 58 | +from nvalchemi.dynamics.hooks import LoggingHook |
| 59 | +from pymatgen.util.testing import PymatgenTest |
| 60 | + |
| 61 | +import matgl |
| 62 | +from matgl.ext.alchmtk import TensorNetWrapper |
| 63 | + |
| 64 | +# %% |
| 65 | +# Load and wrap model |
| 66 | +# -------------------- |
| 67 | +# Load a pretrained TensorNet PES potential from MatGL and wrap it |
| 68 | +# with ``TensorNetWrapper`` for use in nvalchemi dynamics. |
| 69 | + |
| 70 | +potential = matgl.load_model("TensorNet-MatPES-PBE-v2025.1-PES") |
| 71 | +model = TensorNetWrapper.from_potential(potential) |
| 72 | + |
| 73 | +print(f"Model: TensorNet (cutoff={model.model.cutoff} A, units={model.model.units})") |
| 74 | +print(f"Outputs: {model.model_config.outputs}") |
| 75 | +print(f"Active: {model.model_config.active_outputs}") |
| 76 | + |
| 77 | +# %% |
| 78 | +# Build structure |
| 79 | +# ---------------- |
| 80 | +# Load the LiFePO4 structure (28 atoms) from pymatgen's test database. |
| 81 | +# ``from_structure`` handles periodic boundary conditions automatically. |
| 82 | + |
| 83 | +structure = PymatgenTest.get_structure("LiFePO4") |
| 84 | +n_atoms = len(structure) |
| 85 | + |
| 86 | +# %% |
| 87 | +# Build AtomicData with initial fields |
| 88 | +# -------------------------------------- |
| 89 | +# The integrator requires ``forces``, ``energy``, and ``velocities`` |
| 90 | +# to be present on the batch before the first step. We initialize |
| 91 | +# forces and energy to zero (the model overwrites them at the first |
| 92 | +# BEFORE_COMPUTE hook) and sample velocities from Maxwell-Boltzmann. |
| 93 | + |
| 94 | +T_TARGET = 300.0 # K |
| 95 | +KB_EV = 8.617333262e-5 # eV/K |
| 96 | + |
| 97 | +data = AtomicData.from_structure(structure) |
| 98 | +data.forces = torch.zeros(n_atoms, 3) |
| 99 | +data.energy = torch.zeros(1, 1) |
| 100 | + |
| 101 | +# Maxwell-Boltzmann velocities at T_TARGET |
| 102 | +torch.manual_seed(42) |
| 103 | +masses = data.atomic_masses # amu |
| 104 | +v_scale = (KB_EV * T_TARGET / masses).sqrt().unsqueeze(-1) |
| 105 | +velocities = torch.randn(n_atoms, 3) * v_scale |
| 106 | +velocities -= velocities.mean(dim=0, keepdim=True) # zero COM velocity |
| 107 | +data.add_node_property("velocities", velocities) |
| 108 | + |
| 109 | +batch = Batch.from_data_list([data]) |
| 110 | +print(f"\nStructure: {structure.formula} ({n_atoms} atoms, cubic {structure.lattice.a:.1f} A)") |
| 111 | + |
| 112 | +# %% |
| 113 | +# NVTLangevin integrator and hooks |
| 114 | +# ---------------------------------- |
| 115 | +# :class:`~nvalchemi.dynamics.NVTLangevin` implements the BAOAB Langevin |
| 116 | +# splitting scheme. Key arguments: |
| 117 | +# |
| 118 | +# * ``dt`` — timestep in fs |
| 119 | +# * ``temperature`` — target temperature in K |
| 120 | +# * ``friction`` — Langevin friction in 1/fs |
| 121 | +# * ``random_seed`` — reproducible stochastic forces |
| 122 | +# |
| 123 | +# The neighbor list hook is registered via ``model.make_neighbor_hooks()``, |
| 124 | +# which reads the model's ``NeighborConfig`` and creates the appropriate hook. |
| 125 | + |
| 126 | +nvt = NVTLangevin( |
| 127 | + model=model, |
| 128 | + dt=1.0, |
| 129 | + temperature=T_TARGET, |
| 130 | + friction=0.01, |
| 131 | + random_seed=42, |
| 132 | + n_steps=1000, |
| 133 | +) |
| 134 | + |
| 135 | +for hook in model.make_neighbor_hooks(): |
| 136 | + nvt.register_hook(hook, stage=DynamicsStage.BEFORE_COMPUTE) |
| 137 | + |
| 138 | +with LoggingHook(backend="csv", log_path="md_log.csv", frequency=10) as log_hook: |
| 139 | + nvt.register_hook(log_hook) |
| 140 | + |
| 141 | + print(f"\nRunning {nvt.n_steps} NVT steps at T={T_TARGET} K ...") |
| 142 | + batch = nvt.run(batch) |
| 143 | + print(f"NVT completed {nvt.step_count} steps.") |
| 144 | + |
| 145 | +# %% |
| 146 | +# Inspecting temperature |
| 147 | +# ----------------------- |
| 148 | +# The instantaneous kinetic temperature is computed from the |
| 149 | +# equipartition theorem: |
| 150 | +# |
| 151 | +# T = (2 · KE) / (3 · N · kB) |
| 152 | + |
| 153 | +masses = batch.atomic_masses # (N,) amu |
| 154 | +vels = batch.velocities # (N, 3) |
| 155 | +ke_ev = 0.5 * (masses * (vels**2).sum(dim=-1)).sum().item() |
| 156 | +T_final = (2.0 * ke_ev) / (3.0 * n_atoms * KB_EV) |
| 157 | + |
| 158 | +print(f"\nFinal temperature: {T_final:.1f} K (target: {T_TARGET} K)") |
| 159 | +print(f"Final energy: {batch.energy.item():.4f} eV") |
| 160 | +print("MD trajectory saved to md_log.csv") |
0 commit comments