6.5 Program 6.1, Fortran Code
103
real fxi, fyi, fzi, rij, rrijsq
real rijsq ,rxij, ryij, rzij, vij, wij, oij, rpij
c ***************************************************************
c ** zero forces and potential **
c ***************************************************************
do 10 i = 1, n
fx(i) = 0.0
fy(i) = 0.0
fz(i) = 0.0
10 continue
v = 0.0
r = 0.0
c ** loop over all pairs of atoms **
do 100 i = 1, n - 1
rxi = rx(i)
ryi = ry(i)
rzi = rz(i)
fxi = fx(i)
fyi = fy(i)
fzi = fz(i)
ic = 3 * ( i - 1) + 1
do 99 j = i + 1, n
rxij = rxi - rx(j)
ryij = ryi - ry(j)
rzij = rzi - rz(j)
rijsq = rxij * rxij + ryij * ryij + rzij * rzij
c ** calculate off-diagonal blocks of diffusion tensor **
c ** here we assume the rotne-prager tensor form **
c ** take rpij = 0 instead below for oseen tensor **
jc = ( j - 1 ) * 3 + 1
rpij = 0.
rij = sqrt( rijsq )
rrijsq = 1.0 / rijsq
oij = consij / rij
d( ic , jc ) = oij + rpij +
* ( oij - 3.0 * rpij ) * rxij * rxij * rrijsq
d( ic+1, jc+1 ) = oij + rpij +
* ( oij - 3.0 * rpij ) * ryij * ryij * rrijsq
d( ic+2, jc+2 ) = oij + rpij +
* ( oij - 3.0 * rpij ) * rzij * rzij * rrijsq
d( ic , jc+1 ) =
* ( oij - 3.0 * rpij ) * rxij * ryij * rrijsq
d( ic , jc+2 ) =
* ( oij - 3.0 * rpij ) * rxij * rzij * rrijsq
d( ic+1, jc+2 ) =