double ugb (void)
{
int i,j;
double v,rxi,ryi,rzi,uxi,uyi,uzi,v6,v12,rij,v2,vp;
double rxij,ryij,rzij,uij,rui,ruj,rup,rum,hiu,hisu,sigu,epsu,epss;
v=0;
for (i=1;i<=(n-1);i++)
{
rxi=rx[i]; ryi=ry[i]; rzi=rz[i]; uxi=ux[i]; uyi=uy[i]; uzi=uz[i];
for (j=i+1;j<=n;j++)
{
rxij=rxi-rx[j];
if (rxij<=-l2) rxij+=l;
if (rxij>=l2) rxij-=l;
ryij=ryi-ry[j];
if (ryij<=-l2) ryij+=l;
if (ryij>=l2) ryij-=l;
rzij=rzi-rz[j];
if (rzij<=-l2) rzij+=l;
if (rzij>=l2) rzij-=l;
rij=sqrt(rxij*rxij+ryij*ryij+rzij*rzij);
if (rij<rc)
{
if ((i<=n1) && (j<=n1)){hi=hi1;his=his1;}
else
{ if ((i>n1) && (j>n1)) {hi=hi2;his=his2;}
else {hi=hi12;his=his12;}
}
uij=uxi*ux[j]+uyi*uy[j]+uzi*uz[j];
rui=(rxij*uxi+ryij*uyi+rzij*uzi)/rij;
ruj=(rxij*ux[j]+ryij*uy[j]+rzij*uz[j])/rij;
rup=(rui+ruj)*(rui+ruj); rum=(rui-ruj)*(rui-ruj);
hiu=hi*uij; hisu=his*uij;
sigu=1/sqrt(1-(hi/2)*(rup/(1+hiu)+rum/(1-hiu)));
epsu=1/sqrt(1-hi*hi*uij*uij);
epss=1-(his/2)*(rup/(1+hisu)+rum/(1-hisu));
v2=1/(rij-sigu+1);
v2*=v2; v6=v2*v2*v2; v12=v6*v6;
vp=epsu*epss*epss*(v12-v6);
v+=vp;
}
}
}
return(4*v);
}