"""Calculate the vibrational modes of a H2O molecule."""
from Dacapo import Dacapo
from ASE.Utilities.VibrationalCalculation import *

#Change filename to *.nc output file of your relaxation
loa=Dacapo.ReadAtoms('H2Orelax.nc')
loa.GetCalculator().SetNetCDFFile('temp.nc')

vib=VibrationalCalculation(filetoken='H2O_vib',
                           atoms=loa,
                           freeatoms=[0,1,2],
                           displacements=[0.08] * 3,
                           method=0)


vib.RunCalculations()

#Ensure symmetry of Hessian matrix
vib.SetHessianSymmetry(1)

vib.PrintFrequencies()

#Make trajectory files to visualize normal modes
#for all 9 modes
for mode in range(0,9):
    vib.CreateModeTrajectory(mode=mode,scaling=2)


print 'Zero-point energy = %1.2f eV' % vib.GetZeroPointEnergy()