1#!/usr/bin/env python
2
3from Siesta.Interface import  Atom, Crystal
4from Siesta.calculator import SiestaCalculator
5from ASE.Dynamics.MDMin import MDMin
6
7import urllib
8import os, shutil
9
10URLbase = "http://fisica.ehu.es/ag/siesta-psffiles/"
11
12atoms = Crystal([Atom('H', (0.757,0.586,0.),label="H_test"),
13                     Atom('H', (-0.757,0.586,0.),magmom=1),
14                     Atom('O', (0, 0, 0),magmom=1)],
15                     cell=(6.0, 6.0, 6.0), periodic=1)
16
17# create a work subdirectory
18orig_dir = os.getcwd()
19dir = "ase_relax_work"
20if os.path.isdir(dir): # does dir exist?
21  shutil.rmtree(dir) # yes, remove old directory
22os.mkdir(dir) # make dir directory
23os.chdir(dir) # move to dir
24
25urllib.urlretrieve(URLbase + "H.psf", "H_test.psf")
26urllib.urlretrieve(URLbase + "H.psf", "H.psf")
27urllib.urlretrieve(URLbase + "O.psf", "O.psf")
28
29b = SiestaCalculator(executable="/Users/ag/bin/siesta-xlf")
30b.SetOption("DM.Tolerance","0.001")
31atoms.SetCalculator(b)
32
33energy = atoms.GetPotentialEnergy()
34forces = atoms.GetCartesianForces()
35stress = atoms.GetStress()
36
37print "The energy is: ", energy
38print "Forces ", forces
39print "Stress ", stress
40
41dyn = MDMin(atoms, dt=0.08, fmax=0.02)
42dyn.Converge()
43
44energy = atoms.GetPotentialEnergy()
45forces = atoms.GetCartesianForces()
46stress = atoms.GetStress()
47
48print "The energy is: ", energy
49print "Forces ", forces
50print "Stress ", stress
51
52b.stop()
53
54os.chdir(orig_dir)
55