ASE NVE MD simulation of Organic Molecule

ASE_NVE.py

"""Demonstrates molecular dynamics with constant energy."""

from ase.lattice.cubic import FaceCenteredCubic
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase.calculators.dftb import Dftb
from ase.optimize import QuasiNewton,fire,gpmin,mdmin,LBFGS  
from ase.io import write,read
from ase.build import molecule
from rdkit import Chem
from rdkit.Chem import AllChem
from ase.optimize.minimahopping import MinimaHopping
from ase import *
from ase.optimize.basin import BasinHopping
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase import units

# Set the momenta corresponding to T=300K
atoms = read('trail.mol')
atoms.set_calculator(Dftb(label='h2o',
                         atoms=atoms,
                         Hamiltonian_MaxAngularMomentum_='',
                         Hamiltonian_MaxAngularMomentum_C='"p"',
                         Hamiltonian_MaxAngularMomentum_N='"p"',
                         Hamiltonian_MaxAngularMomentum_O='"p"',
                         Hamiltonian_MaxAngularMomentum_H='"s"',
                         ))
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)

# We want to run MD with constant energy using the VelocityVerlet algorithm.
dyn = VelocityVerlet(atoms, 0.5 * units.fs)  # 5 fs time step.


def printenergy(a=atoms):  # store a reference to atoms in the definition.
    """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))

# Now run the dynamics
dyn.attach(printenergy, interval=10)
printenergy()
dyn.run(200)