"docs/vscode:/vscode.git/clone" did not exist on "969e2af866045417dccbc3980422c80d9736d970"
ase_interface.py 2.31 KB
Newer Older
Gao, Xiang's avatar
Gao, Xiang committed
1
2
# -*- coding: utf-8 -*-
"""
3
4
Structure minimization and constant temperature MD using ASE interface
======================================================================
Gao, Xiang's avatar
Gao, Xiang committed
5

6
7
8
This example is modified from the official `home page` and
`Constant temperature MD`_ to use the ASE interface of TorchANI as energy
calculator.
Gao, Xiang's avatar
Gao, Xiang committed
9

10
11
.. _home page:
    https://wiki.fysik.dtu.dk/ase/
Gao, Xiang's avatar
Gao, Xiang committed
12
13
14
15
16
17
18
19
20
.. _Constant temperature MD:
    https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html#constant-temperature-md
"""


###############################################################################
# To begin with, let's first import the modules we will use:
from ase.lattice.cubic import Diamond
from ase.md.langevin import Langevin
21
from ase.optimize import BFGS
Gao, Xiang's avatar
Gao, Xiang committed
22
23
24
25
26
27
28
from ase import units
import torchani


###############################################################################
# Now let's set up a crystal
atoms = Diamond(symbol="C", pbc=True)
29
print(len(atoms), "atoms in the cell")
Gao, Xiang's avatar
Gao, Xiang committed
30
31
32

###############################################################################
# Now let's create a calculator from builtin models:
Gao, Xiang's avatar
Gao, Xiang committed
33
calculator = torchani.models.ANI1ccx().ase()
Gao, Xiang's avatar
Gao, Xiang committed
34
35
36
atoms.set_calculator(calculator)

###############################################################################
37
38
39
40
41
# Now let's minimize the structure:
print("Begin minimizing...")
opt = BFGS(atoms)
opt.run(fmax=0.001)
print()
Gao, Xiang's avatar
Gao, Xiang committed
42
43
44


###############################################################################
45
46
# Now create a callback function that print interesting physical quantities:
def printenergy(a=atoms):
Gao, Xiang's avatar
Gao, Xiang committed
47
48
49
50
51
52
53
    """Function to print the potential, kinetic and total energy."""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))


54
55
56
57
58
###############################################################################
# We want to run MD with constant energy using the Langevin algorithm
# with a time step of 1 fs, the temperature 300K and the friction
# coefficient to 0.02 atomic units.
dyn = Langevin(atoms, 1 * units.fs, 300 * units.kB, 0.2)
59
dyn.attach(printenergy, interval=50)
Gao, Xiang's avatar
Gao, Xiang committed
60
61
62

###############################################################################
# Now run the dynamics:
63
print("Beginning dynamics...")
Gao, Xiang's avatar
Gao, Xiang committed
64
65
printenergy()
dyn.run(500)