Add minimization, and choose better parameter for dynamics for the ASE example (#127)

examples/ase_langevin.py → examples/ase_interface.py

 # -*- coding: utf-8 -*-
 """
-Constant temperature MD using ASE interface
-===========================================
+Structure minimization and constant temperature MD using ASE interface
+======================================================================
 
-This example is modified from the official `Constant temperature MD`_ to use
-the ASE interface of TorchANI as energy calculator.
+This example is modified from the official `home page` and
+`Constant temperature MD`_ to use the ASE interface of TorchANI as energy
+calculator.
 
+.. _home page:
+    https://wiki.fysik.dtu.dk/ase/
 .. _Constant temperature MD:
     https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html#constant-temperature-md
 """
@@ -15,6 +18,7 @@ the ASE interface of TorchANI as energy calculator.
 # To begin with, let's first import the modules we will use:
 from ase.lattice.cubic import Diamond
 from ase.md.langevin import Langevin
+from ase.optimize import BFGS
 from ase import units
 import torchani
 
@@ -22,6 +26,7 @@ import torchani
 ###############################################################################
 # Now let's set up a crystal
 atoms = Diamond(symbol="C", pbc=True)
+print(len(atoms), "atoms in the cell")
 
 ###############################################################################
 # Now let's create a calculator from builtin models:
@@ -31,15 +36,16 @@ calculator = torchani.ase.Calculator(builtin.species, builtin.aev_computer,
 atoms.set_calculator(calculator)
 
 ###############################################################################
-# We want to run MD with constant energy using the Langevin algorithm
-# with a time step of 5 fs, the temperature 50K and the friction
-# coefficient to 0.02 atomic units.
-dyn = Langevin(atoms, 5 * units.fs, 50 * units.kB, 0.002)
+# Now let's minimize the structure:
+print("Begin minimizing...")
+opt = BFGS(atoms)
+opt.run(fmax=0.001)
+print()
 
 
 ###############################################################################
-# Let's print energies every 50 steps:
-def printenergy(a=atoms):  # store a reference to atoms in the definition.
+# Now create a callback function that print interesting physical quantities:
+def printenergy(a=atoms):
     """Function to print the potential, kinetic and total energy."""
     epot = a.get_potential_energy() / len(a)
     ekin = a.get_kinetic_energy() / len(a)
@@ -47,9 +53,15 @@ def printenergy(a=atoms):  # store a reference to atoms in the definition.
           'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))
 
 
-dyn.attach(printenergy, interval=50)
+###############################################################################
+# 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)
+dyn.attach(printenergy, interval=500)
 
 ###############################################################################
 # Now run the dynamics:
+print("Beginning dynamics...")
 printenergy()
 dyn.run(500)