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