[ase-users] emt giving nonzero interaction energies at large separation

[Tom Haxton]

Hello,
I want to use ase to calculate the interaction energy between a metal slab
and a substituted benzene molecule as a function of relative position.  I
expect that the interaction energy should go to zero as the separation
between the molecules increases.  However, using the emt calculator, this is
not the case.  At large separations, the sum of the energies of the isolated
slab and the isolated molecule is not equal to the energy of the combined
system.  I use the same cell size, shape, and periodic boundary conditions
in all three calculations.  For large separations and large cell sizes, the
result is independent of separation, cell size, and boundary conditions.

I am using ase-3.4.  I also tried ase-3.5 and got the same problem, though a
different number for the energy discrepancy.  I am working on a x86_64
machine running Linux 2.6, python 2.6.5, and numpy 1.3.0.

Below is the script that produces the error.  It is the more minimal case of
a single isolated Au atom interacting with a Benzene molecule, separated by
50 angstroms, in a 100 x 100 x 100 angstrom cell.  The energies are:

Benzene    -3.47792756546
Au                3.8
Combined    -1.46565762747
Difference   -1.78773006201

Can anyone tell me if there is something I am doing wrong or if this is a
bug?  Also, is there any way to get out any useful intermediate information
from the emt calculation?  Does ASAP have the same problem?  (I don't have
ASAP installed yet.)  Thank you very much.

Best,
Tom



#/usr/bin/python

import sys
import math
import copy
from ase import Atoms
from ase.visualize import view
from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton
from ase.lattice.surface import fcc111,add_adsorbate

mycalculator = EMT()
mybigcell = (100., 100., 100.)
separation = 50.

r_CC = 1.374
r_HC = 1.084

benzene = Atoms('CCCCCCHHHHHH', positions=[(r_CC, 0., 0.), (-r_CC, 0., 0),
(0.5 * r_CC, 0.5 * math.sqrt(3.0) * r_CC, 0.), (0.5 * r_CC, -0.5 *
math.sqrt(3.0) * r_CC, 0.), (-0.5 * r_CC, 0.5 * math.sqrt(3.0) * r_CC, 0.),
(-0.5 * r_CC, -0.5 * math.sqrt(3.0) * r_CC, 0.), (r_CC + r_HC, 0., 0.),
(-(r_CC + r_HC), 0., 0.), (0.5 * (r_CC + r_HC), 0.5 * math.sqrt(3.0) * (r_CC
+ r_HC), 0), (0.5 * (r_CC + r_HC), -0.5 * math.sqrt(3.0) * (r_CC + r_HC),
0), (-0.5 * (r_CC + r_HC), 0.5 * math.sqrt(3.0) * (r_CC + r_HC), 0), (-0.5 *
(r_CC + r_HC), -0.5 * math.sqrt(3.0) * (r_CC + r_HC), 0)])

gold = Atoms('Au', positions=[(0., 0., 0.)])

molecule1 = copy.deepcopy(benzene)
molecule2 = copy.deepcopy(gold)

molecule1.set_cell(mybigcell)
molecule1.set_pbc(False)
molecule1.set_calculator(mycalculator)
molecule2.set_cell(mybigcell)
molecule2.set_pbc(False)
molecule2.set_calculator(mycalculator)
molecule2.translate((0., 0., separation))
moleculecombined = molecule1 + molecule2
molecule1.center()
molecule2.center()
moleculecombined.center()
view(molecule1)
view(molecule2)

e1 = molecule1.get_potential_energy()
e2 = molecule2.get_potential_energy()

moleculecombined.set_cell(mybigcell)
moleculecombined.set_pbc(False)
moleculecombined.set_calculator(mycalculator)

ecombined = moleculecombined.get_potential_energy()

print ecombined - e1 - e2, e1, e2, ecombined

#view(moleculecombined)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20110912/38e97ce8/attachment.html>