C
C $Id: pme_calcb.F,v 1.2 1998/07/16 16:40:13 jjv5 Exp $
C
C------------------------------------------------------------------------
      subroutine pme_calcb

      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     this routine calculates |1/bi|**2, i=x,y,z as given by [Eq.4.4]
C     of Essmann et.al., J.Chem.Phys. 103 8577
C
C     written by Arjan van der Vaart

      dimension smn(maxspline)

CC-------------------------------------------------------------------CC

C     some constants
      pi2k3pme = pi2/dk3pme
      pi2k2pme = pi2/dk2pme
      pi2k1pme = pi2/dk1pme

C     do the splining
      call mkbspline(smn)

C     |bz|**2
      do 20 m3=0,k3pmem1
         argz = pi2k3pme*dble(m3)
         bzr = smn(1)
         bzi = 0.0
         do 10 k=1,nspline-2
            arg = argz*dble(k)
            bzr = bzr + smn(k+1)*cos(arg)
            bzi = bzi + smn(k+1)*sin(arg)
 10      continue
         bfaczpme(m3) = bzr*bzr + bzi*bzi
 20   continue

C     |by|**2
      do 40 m2=0,k2pmem1
         argy = pi2k2pme*dble(m2)
         byr = smn(1)
         byi = 0.0
         do 30 k=1,nspline-2
            arg = argy*dble(k)
            byr = byr + smn(k+1)*cos(arg)
            byi = byi + smn(k+1)*sin(arg)
 30      continue
         bfacypme(m2) = byr*byr + byi*byi
 40   continue

C     |bx|**2
      do 60 m1=0,k1pmem1
         argx = pi2k1pme*dble(m1)
         bxr = smn(1)
         bxi = 0.0
         do 50 k=1,nspline-2
            arg = argx*dble(k)
            bxr = bxr + smn(k+1)*cos(arg)
            bxi = bxi + smn(k+1)*sin(arg)
 50      continue
         bfacxpme(m1) = bxr*bxr + bxi*bxi
 60   continue

      end