double press (void)
{
int i,j;
double p,rxi,ryi,rzi,fij,uxi,uyi,uzi;
double rxij,ryij,rzij,uij,rui,ruj,rup,rum,hiu,hisu,sigu,epsu,epss;
double rij,v1,v2,v6,v7,v12,v13,pvirial;
double dudx,dudy,dudz,drrdx,drrdy,drrdz,dgdx,dgdy,dgdz,dgsdx,dgsdy,dgsdz;
p=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));
v1=1/(rij-sigu+1);
v2=v1*v1;
v6=v2*v2*v2;
v7=v6*v1;
v12=v6*v6;
v13=v12*v1;
dgdx=-hi/rij*((rui+ruj)/(1+hiu)*(uxi+ux[j])+(rui-ruj)/(1-hiu)*(uxi-ux[j]));
dgdx+=rxij*hi/(rij*rij)*(rup/(1+hiu)+rum/(1-hiu));
dgsdx=(rui+ruj)/(1+hisu)*(uxi+ux[j])+(rui-ruj)/(1-hisu)*(uxi-ux[j]);
dgsdx*=-his/rij;
dgsdx+=rxij*his/(rij*rij)*(rup/(1+hisu)+rum/(1-hisu));
drrdx=rxij/rij+0.5*sigu*sigu*sigu*dgdx;
dudx=4*epsu*epss*epss*(-12*v13+6*v7)*drrdx+4*(v12-v6)*2*epss*dgsdx;
dgdy=-hi/rij*((rui+ruj)/(1+hiu)*(uyi+uy[j])+(rui-ruj)/(1-hiu)*(uyi-uy[j]));
dgdy+=ryij*hi/(rij*rij)*(rup/(1+hiu)+rum/(1-hiu));
dgsdy=(rui+ruj)/(1+hisu)*(uyi+uy[j])+(rui-ruj)/(1-hisu)*(uyi-uy[j]);
dgsdy*=-his/rij;
dgsdy+=ryij*his/(rij*rij)*(rup/(1+hisu)+rum/(1-hisu));
drrdy=ryij/rij+0.5*sigu*sigu*sigu*dgdy;
dudy=4*epsu*epss*epss*(-12*v13+6*v7)*drrdy+4*(v12-v6)*2*epss*dgsdy;
dgdz=-hi/rij*((rui+ruj)/(1+hiu)*(uzi+uz[j])+(rui-ruj)/(1-hiu)*(uzi-uz[j]));
dgdz+=rzij*hi/(rij*rij)*(rup/(1+hiu)+rum/(1-hiu));
dgsdz=(rui+ruj)/(1+hisu)*(uzi+uz[j])+(rui-ruj)/(1-hisu)*(uzi-uz[j]);
dgsdz*=-his/rij;
dgsdz+=rzij*his/(rij*rij)*(rup/(1+hisu)+rum/(1-hisu));
drrdz=rzij/rij+0.5*sigu*sigu*sigu*dgdz;
dudz=4*epsu*epss*epss*(-12*v13+6*v7)*drrdz+4*(v12-v6)*2*epss*dgsdz;
pvirial=-rxij*dudx-ryij*dudy-rzij*dudz;
p+=pvirial;
}
}
}
return (ro*t+p/(3*l*l*l));
}