This is the mail archive of the
gcc-bugs@gcc.gnu.org
mailing list for the GCC project.
This code runs *very* slowly
- To: egcs-bugs at cygnus dot com
- Subject: This code runs *very* slowly
- From: Rahul Siddharthan <rsidd at physics dot iisc dot ernet dot in>
- Date: Mon, 2 Nov 1998 13:07:20 +0530 (IST)
Hi
I'm reporting a problem of extreme slowness of the egcs 1.1b C
compiler, when running some code I have which has rather unusual
requirements. This particular code runs four times slower with
egcs than with Digital's C compiler, and also takes several
minutes to compile, while Digital's compiler compiles it in
seconds. Most other code runs at a comparable speed (perhaps 20%
slower) to code produced by Digital's compiler.
I'm a physics student running a monte-carlo simulation code for
some spin systems. I have to model some rather complicated
nearest-neighbour, next-nearest neighbour, etc, interactions, and
while the nearest neighbour isn't too hard, the further ones are
very numerous. I decided the most efficient way to do it was to
write out each interaction other than the near-neighbour ones
explicitly, define them separately in a header file, and
"include" it whenever I want to evaluate the energy. The code is
below, but it will be a bit hard to explain without going into
the physics of the problem. Probably a lot of it could be better
written but I'm not a very experienced programmer.
We have four Digital Alpha workstations, two of which (333 Mhz)
run Digital Unix 4.0 with Digital's c compiler, and two (433 MHz)
run Red Hat Linux 5.0 with egcs 1.1b. I tried compiling on both,
just to compare the speed. When I have only near-neighbour
interactions, I found that the Digital C code runs around 10%-20%
faster than the egcs code -- despite the machine clock speed
being around 25% less. This speed difference was acceptable to me
since I dont expect gcc to be as fast as a compiler specifically
optimized for the Alpha architecture.
When I tried including the next-neighbour interactions through
#defines, the Digital C code slowed a bit (as expected), but the
egcs code slowed an enormous lot: it now runs around a factor of
4 slower than the Digital code. This is way unacceptable: even
if there are three other jobs running on the Digital Unix machine
it still runs faster than on an otherwise idle Linux machine.
Including more further interactions made the performance still
worse. The compile time also grows greatly: with addition of
beyond-second-neighbour interactions the code takes over half an
hour to compile, while Digital's compiler takes less than a
minute for all three cases.
The command line I used was
(Digital)
cc -O5 -lm -D_BSD -I. pyromc.c
(egcs)
gcc -O -lm -I. pyromc.c
With egcs I also experimented with -O2 (which in fact makes it
slightly faster), -funroll-loops and -ffast-math (which have no
particular difference).
The actual code follows. This is in the next neighbour case. The
only difference for nearest-neighbours is, remove the #defines
for dhseq1, dhseq2, dhseq3, dhseq4 in energyseq.h and all
references to them and the function energylatticenn in the main
program. I presumed that #defining each of the next-neighbour
interactions as a macro, rather than calling a function to
evaluate each separately (I could need between 24 and 300 of
each) would be much faster, and that indeed seems the case with
Digital's compiler, but with gcc there seems to be some other
overhead.
Rahul Siddharthan rsidd@physics.iisc.ernet.in
----------------------------------------------------------------------------
Department of Physics, Indian Institute of Science, Bangalore 560012 (India)
Phone: (Off) (91-80)3092596 (Hostel) 3092624, room U-70 Fax:(91-80)3341683
pyromc.c:
--------------------------
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <energyseq.h>
#define jsign 1
#define nn 5
#define lmax 10
#define randomint(a) floor(random()/(RAND_MAX+1.0)*a)
main()
{
short int a,b,c,lx,ly,lz;
unsigned int arrindex;
unsigned long int ctr,n;
unsigned int energies[2400];
double E,T,Eavg,Etot;
double energycell[16];
short int lattice[lmax][lmax][lmax];
short int randstate[256];
unsigned int randseed;
int randn;
double energylatticen();
double energylatticenn();
double newstate();
FILE *ptvar;
randseed=134;
randn=256;
initstate(randseed,randstate,randn);
makeenergycell(energycell);
for (a=0; a<10; a++) {
for (b=0; b<10; b++) {
for (c=0; c<10; c++) {
lattice[a][b][c]=12;
}
}
}
lx=4;
ly=4;
lz=4;
E=energylatticen(lx,ly,lz,lattice,energycell)+energylatticenn(lx,ly,lz,lattice);
printf("Energy is %f\n",E);
for (T=3.0; T>=0.7; T=T-0.005){
Etot=0;
for (n=1; n<=1000; n++) {
E=newstate(lx,ly,lz,lattice,energycell,E,T,1000,&Eavg);
if (n>400) Etot=Etot+Eavg;
}
Etot=Etot/600;
ptvar=fopen("datafile","a");
fprintf(ptvar,"%f %f\n",T,Etot);
fclose(ptvar);
}
}
/*
This simulates a pyrochlore lattice which is generated by
translating a tetrahedron with vertices at
(consult notes)
through lattice translation vectors
()
ie a lattice with four basis atoms.
Each unit cell then consists of a tetrahedron, the vertices of
which have Ising spins pointing out or into the tetrahedron.
Represent outward spin by 1, inward by 0. Then the configuration of
a tetrahedron can be given by a four-bit integer.
a1 a2 a3 a4
=
a3
a4
a1 a2
where a4 is the apex atom (above the plane of a1,a2,a3).
energy1 represents interaction of the tetrahedron with the next
tetrahedron given by lattice translation vector 1,
energy2, energy3 likewise.
*/
makeenergycell(energycell)
/* This is the energy of one tetrahedron. */
double energycell[];
{
short int a;
for (a=0; a<=15; a++) {
if ((a==3)||(a==5)||(a==6)||(a==9)||(a==10)||(a==12))
energycell[a]=-2*jsign;
else if ((a==0) || (a==15))
energycell[a]=6*jsign;
else
energycell[a]=0*jsign;
}
}
double energylatticen(lx,ly,lz,lattice,energycell)
short int lx,ly,lz;
short int lattice[][lmax][lmax];
double energycell[];
{
double elat;
short int l1,l2,l3,l2m,l3m,a,a1,a2,a3,a4,a5,a6;
elat=0;
for (l1=0; l1<lx; l1++) {
for (l2=0; l2<ly; l2++) {
for (l3=0; l3<lz; l3++) {
l2m=(l2-1+lx)%lx;
l3m=(l3-1+lx)%lx;
a=lattice[l1][l2][l3];
a1=lattice[(l1+1)%lx][l2][l3];
a2=lattice[l1][(l2+1)%ly][l3];
a3=lattice[l1][l2][(l3+1)%lz];
a4=lattice[(l1+1)%lx][l2m][l3];
a5=lattice[(l1+1)%lx][l2][l3m];
a6=lattice[l1][(l2+1)%ly][l3m];
elat=elat+energycell[a]+energy1(a,a1)+energy2(a,a2)+energy3(a,a3)
+energy4(a,a4)+energy5(a,a5)+energy6(a,a6);
}
}
}
return elat;
}
double energylatticenn(lx,ly,lz,lattice)
short int lx,ly,lz;
short int lattice[][lmax][lmax];
{
double elat,dh;
short int l1,l2,l3,l2m,l3m,a,a1,a2,a3,a4,a5,a6,an1;
elat=0;
for (l1=0; l1<lx; l1++) {
for (l2=0; l2<ly; l2++) {
for (l3=0; l3<lz; l3++) {
dh=0.0;
a=lattice[l1][l2][l3];
dhseq1;
dhseq2;
dhseq3;
dhseq4;
elat=elat+dh/2;
}
}
}
elat=elat*nn;
return elat;
}
double newstate(lx,ly,lz,lattice,energycell,E,T,Niters,Eavg)
short int lx,ly,lz;
double T,E,*Eavg;
double energycell[];
short int lattice[][lmax][lmax];
unsigned long int Niters;
/* Generates a new configuration by randomly changing the configuration
of a randomly chosen tetrahedron, and does this Niter times */
{
short int anew,a,aold,atmp,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,l1,l2,l3,
n,randint,an1;
unsigned long int m;
double dhn,dh,en1,en2,entot,tw,twrand;
en1=E; entot=E;
for (m=0; m<Niters; m++) {
l1=randomint(lx);
l2=randomint(ly);
l3=randomint(lz); /* Coordinates of randomly chosen tetrahedron */
a= lattice[l1][l2][l3];
aold=a;
randint=randomint(4);
atmp=1 << randint; /* New configuration of this tetrahedron */
anew = (a&~atmp) | (~a&atmp);
if (atmp==8){
a7=lattice[(l1-1+lx)%lx][l2][l3];
a10=lattice[(l1-1+lx)%lx][(l2+1)%ly][l3];
a11=lattice[(l1+lx-1)%lx][l2][(l3+1)%lz];
dhn=energycell[a] +energy1(a7,a) +energy4(a10,a)+energy5(a11,a);
dh=0;
dhseq1;
a=anew;
dhn=-dhn+energycell[a] +energy1(a7,a) +energy4(a10,a)+energy5(a11,a);
dh=-dh;
dhseq1;
}
else
if (atmp==4){
a4=lattice[(l1+1)%lx][(l2-1+ly)%ly][l3];
a8=lattice[l1][(l2-1+ly)%ly][l3];
a12=lattice[l1][(l2-1+ly)%ly][(l3+1)%lz];
dhn=energycell[a] +energy4(a,a4) +energy2(a8,a) +energy6(a12,a);
dh=0;
dhseq2;
a=anew;
dhn=-dhn+energycell[a] +energy4(a,a4) +energy2(a8,a) +energy6(a12,a);
dh=-dh;
dhseq2;
}
else
if (atmp==2){
a1=lattice[(l1+1)%lx][l2][l3];
a2=lattice[l1][(l2+1)%ly][l3];
a3=lattice[l1][l2][(l3+1)%lz];
dhn=energycell[a]+energy1(a,a1)+energy2(a,a2)+energy3(a,a3);
dh=0;
dhseq3;
a=anew;
dhn=-dhn+ energycell[a]+energy1(a,a1)+energy2(a,a2)+energy3(a,a3);
dh=-dh;
dhseq3;
}
else
{
a5=lattice[(l1+1)%lx][l2][(l3-1+lz)%lz];
a6=lattice[l1][(l2+1)%ly][(l3-1+lz)%lz];
a9=lattice[l1][l2][(l3-1+lz)%lz];
dhn=energycell[a]+ energy5(a,a5)+energy6(a,a6)+energy3(a9,a);
dh=0;
dhseq4;
a=anew;
dhn= -dhn+energycell[a]+ energy5(a,a5)+energy6(a,a6)+energy3(a9,a);
dh=-dh;
dhseq4;
}
dh=dh*nn;
if ((dh+dhn)>0) {
tw=exp(-(dh+dhn)/T);
twrand=random()/(RAND_MAX+1.0);
if (twrand<=tw) {
en1=en1+dh+dhn;
lattice[l1][l2][l3]=anew;
}
}
else {
en1=en1+dh+dhn;
lattice[l1][l2][l3]=anew;
}
entot=entot+en1;
}
*Eavg=entot/(Niters+1.0);
/*
anew=lattice[l1][l2][l3];
if (anew==aold) {dh=0; dhn=0;}
lattice[l1][l2][l3]=aold;
en1=energylatticen(lx,ly,lz,lattice,energycell)
+energylatticenn(lx,ly,lz,lattice);
lattice[l1][l2][l3]=anew;
en2=energylatticen(lx,ly,lz,lattice,energycell)
+energylatticenn(lx,ly,lz,lattice);
if (floor(en1+dh+dhn) != floor(en2))
printf("The energy now is %f (updated) and %f (exact) \n",en1+dh+dhn,en2);
*/
return en1;
}
--------------------------------------
energyseq.h for next-to-nearest-neighbour interactions
#define energy1(a,b) jsign*(1-(((a>>1) + (b>>3))&1)*2) /* 3,1 */
#define energy2(a,b) jsign*(1-(((a>>1) + (b>>2))&1)*2) /* 3,2 */
#define energy3(a,b) jsign*(1-(((a>>1) + b)&1)*2) /* 3,4 */
#define energy4(a,b) jsign*(1-(((a>>2) + (b>>3))&1)*2) /* 2,1 */
#define energy5(a,b) jsign*(1-((a + (b>>3))&1)*2) /* 4,1 */
#define energy6(a,b) jsign*(1-((a + (b>>2))&1)*2) /* 4,2 */
#define energy7(a,b) jsign*(1-(((a>>3) + (b>>3))&1)*2) /* 1,1 */
#define energy8(a,b) jsign*(1-(((a>>2) + (b>>2))&1)*2) /* 2,2 */
#define energy9(a,b) jsign*(1-(((a>>1) + (b>>1))&1)*2) /* 3,3 */
#define energy10(a,b) jsign*(1-((a + b)&1)*2) /* 4,4 */
#define energy11(a,b) jsign*(1-(((a>>3) + (b>>2))&1)*2) /* 1,2 */
#define energy12(a,b) jsign*(1-(((a>>3) + (b>>1))&1)*2) /* 1,3 */
#define energy13(a,b) jsign*(1-(((a>>3) + b)&1)*2) /* 1,4 */
#define energy14(a,b) jsign*(1-(((a>>2) + (b>>1))&1)*2) /* 2,3 */
#define energy15(a,b) jsign*(1-(((a>>2) + b)&1)*2) /* 2,4 */
#define energy16(a,b) jsign*(1-((a + (b>>1))&1)*2) /* 4,3 */
#define dhseq1 \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy11(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy13(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy11(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy12(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy12(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy13(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy12(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy13(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy12(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy13(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy11(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy11(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*-0.075000;
#define dhseq2 \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy14(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy15(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy4(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy15(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy4(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy14(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy14(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy15(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy14(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy15(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy4(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy4(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*0.075000;
#define dhseq3 \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy3(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy2(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy3(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy1(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy2(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy2(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy1(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy3(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy1(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy1(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy2(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy3(a,an1)*-0.115470;
#define dhseq4 \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy16(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz] ; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy6(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy16(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz] ; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy5(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy6(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz] ; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy5(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy16(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy6(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy5(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy6(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy16(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz] ; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy5(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz] ; \
dh=dh+energy10(a,an1)*0.075000;
--------------------------------
energyseq.h for further interactions
#define energy1(a,b) jsign*(1-(((a>>1) + (b>>3))&1)*2) /* 3,1 */
#define energy2(a,b) jsign*(1-(((a>>1) + (b>>2))&1)*2) /* 3,2 */
#define energy3(a,b) jsign*(1-(((a>>1) + b)&1)*2) /* 3,4 */
#define energy4(a,b) jsign*(1-(((a>>2) + (b>>3))&1)*2) /* 2,1 */
#define energy5(a,b) jsign*(1-((a + (b>>3))&1)*2) /* 4,1 */
#define energy6(a,b) jsign*(1-((a + (b>>2))&1)*2) /* 4,2 */
#define energy7(a,b) jsign*(1-(((a>>3) + (b>>3))&1)*2) /* 1,1 */
#define energy8(a,b) jsign*(1-(((a>>2) + (b>>2))&1)*2) /* 2,2 */
#define energy9(a,b) jsign*(1-(((a>>1) + (b>>1))&1)*2) /* 3,3 */
#define energy10(a,b) jsign*(1-((a + b)&1)*2) /* 4,4 */
#define energy11(a,b) jsign*(1-(((a>>3) + (b>>2))&1)*2) /* 1,2 */
#define energy12(a,b) jsign*(1-(((a>>3) + (b>>1))&1)*2) /* 1,3 */
#define energy13(a,b) jsign*(1-(((a>>3) + b)&1)*2) /* 1,4 */
#define energy14(a,b) jsign*(1-(((a>>2) + (b>>1))&1)*2) /* 2,3 */
#define energy15(a,b) jsign*(1-(((a>>2) + b)&1)*2) /* 2,4 */
#define energy16(a,b) jsign*(1-((a + (b>>1))&1)*2) /* 4,3 */
#define dhseq1 \
an1=lattice[(l1+lx-2)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy12(a,an1)*0.016970; \
an1=lattice[(l1+lx-2)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy13(a,an1)*0.016970; \
an1=lattice[(l1+lx-2)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy11(a,an1)*0.016970; \
an1=lattice[(l1+lx-2)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy12(a,an1)*0.016970; \
an1=lattice[(l1+lx-2)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy11(a,an1)*0.016970; \
an1=lattice[(l1+lx-2)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy13(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy12(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*-0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy12(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy13(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+2)%lz]; \
dh=dh+energy13(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy12(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy11(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy13(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy11(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy12(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+2)%lz]; \
dh=dh+energy13(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy11(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy12(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy12(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy13(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy11(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy13(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+2)%ly][(l3+lz-1)%lz]; \
dh=dh+energy11(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+2)%ly][(l3+lz+0)%lz]; \
dh=dh+energy11(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy11(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy12(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy13(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy11(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy12(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy13(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy11(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy12(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy13(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy11(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy12(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy13(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy11(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy12(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy13(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy11(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy12(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy13(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy12(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy13(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy7(a,an1)*0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy11(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy12(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy7(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy11(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy13(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy7(a,an1)*-0.000000;
#define dhseq2 \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy14(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*-0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy14(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy15(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+2)%lz]; \
dh=dh+energy15(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy4(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy14(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy15(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy4(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy14(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy15(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy14(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy15(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*0.000000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-2)%ly][(l3+lz+1)%lz]; \
dh=dh+energy14(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-2)%ly][(l3+lz+1)%lz]; \
dh=dh+energy15(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy14(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy4(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy15(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy4(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy14(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+2)%lz]; \
dh=dh+energy15(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy4(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy14(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy15(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy4(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy14(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy15(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy4(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy14(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy4(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy15(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-2)%ly][(l3+lz+0)%lz]; \
dh=dh+energy4(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-2)%ly][(l3+lz+0)%lz]; \
dh=dh+energy14(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-2)%ly][(l3+lz+1)%lz]; \
dh=dh+energy4(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-2)%ly][(l3+lz+1)%lz]; \
dh=dh+energy15(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy4(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy14(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy14(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy15(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy4(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy8(a,an1)*0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy15(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy4(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy14(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy15(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy4(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy8(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy14(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy15(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy8(a,an1)*-0.000000; \
an1=lattice[(l1+lx+2)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy4(a,an1)*0.016970; \
an1=lattice[(l1+lx+2)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy4(a,an1)*0.016970;
#define dhseq3 \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy2(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy3(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy1(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy2(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy3(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+2)%lz]; \
dh=dh+energy3(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*-0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy1(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy2(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy3(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy2(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*-0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy3(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+2)%ly][(l3+lz+0)%lz]; \
dh=dh+energy2(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy1(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy3(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy1(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy2(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy3(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+2)%lz]; \
dh=dh+energy3(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy1(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy2(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy1(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy2(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy1(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy2(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy3(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy1(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy3(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy2(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy3(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+2)%ly][(l3+lz-1)%lz]; \
dh=dh+energy2(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*-0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy1(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy2(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy3(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy1(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy9(a,an1)*-0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy3(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy1(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy2(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy3(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy2(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy9(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy3(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy1(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy3(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy1(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy2(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy9(a,an1)*0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy1(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy2(a,an1)*0.016970; \
an1=lattice[(l1+lx+2)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy1(a,an1)*0.016970; \
an1=lattice[(l1+lx+2)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy1(a,an1)*0.016970;
#define dhseq4 \
an1=lattice[(l1+lx-1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy10(a,an1)*0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy16(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy5(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy6(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy16(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy6(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy16(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy6(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy16(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy10(a,an1)*0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy5(a,an1)*0.016970; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy6(a,an1)*-0.115470; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy16(a,an1)*0.003578; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy10(a,an1)*0.000000; \
an1=lattice[(l1+lx-1)%lx][(l2+ly+2)%ly][(l3+lz-1)%lz]; \
dh=dh+energy6(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy16(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy5(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy6(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy16(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy5(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy16(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy5(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy6(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy5(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy6(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+0)%ly][(l3+lz+1)%lz]; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-2)%lz]; \
dh=dh+energy6(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-2)%lz]; \
dh=dh+energy16(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy5(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy16(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy5(a,an1)*0.003578; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy6(a,an1)*-0.115470; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy16(a,an1)*0.016970; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx+0)%lx][(l2+ly+2)%ly][(l3+lz-1)%lz]; \
dh=dh+energy6(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy5(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy16(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy10(a,an1)*0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy5(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy6(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy16(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+0)%lz]; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly-1)%ly][(l3+lz+1)%lz]; \
dh=dh+energy10(a,an1)*0.000000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-2)%lz]; \
dh=dh+energy5(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-2)%lz]; \
dh=dh+energy16(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy6(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy16(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy10(a,an1)*-0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy5(a,an1)*-0.115470; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy6(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy16(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+0)%ly][(l3+lz+0)%lz]; \
dh=dh+energy10(a,an1)*0.075000; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-2)%lz]; \
dh=dh+energy5(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-2)%lz]; \
dh=dh+energy6(a,an1)*0.016970; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy5(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy6(a,an1)*0.003578; \
an1=lattice[(l1+lx+1)%lx][(l2+ly+1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy10(a,an1)*0.000000; \
an1=lattice[(l1+lx+2)%lx][(l2+ly-1)%ly][(l3+lz-1)%lz]; \
dh=dh+energy5(a,an1)*0.016970; \
an1=lattice[(l1+lx+2)%lx][(l2+ly+0)%ly][(l3+lz-1)%lz]; \
dh=dh+energy5(a,an1)*0.016970;
--------------------------------------------