      subroutine mespmm(nreal,ninteger,iptsflag,nauxiliary,nunique,
     &                  npoints,ncentr,ndim,nbsystems,nconts,ncontp,
     &                  ncontd,ncds,ncdspd,nshels,nshelp,nsheld,mtloca,
     &                  iprimvc,nprimvc,icentvc,icdcfunc,icdlfunc,
     &                  nptsvec,isys,iwkvec,coord,charge,alpha,coeff,
     &                  alphacd,coeffscd,coeffpcd,coeffdcd,cdfitc,
     &                  density,tdensity,qqa,qqb,qqd,esppt,wkvec)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      dimension g(6)

      dimension nshels(*),nshelp(*),nsheld(*),mtloca(*),iprimvc(*),
     &          nprimvc(*),icentvc(*),icdcfunc(*),icdlfunc(*),
     &          nptsvec(*),isys(ndim,*),iwkvec(*)

      dimension coord(3,*),charge(*),alpha(*),coeff(*),alphacd(*),
     &          coeffscd(*),coeffpcd(*),coeffdcd(*),cdfitc(*),
     &          density(*),tdensity(*),qqa(2*ndim,*),qqb(*),
     &          qqd(npoints,*),esppt(npoints,*),wkvec(*)

      data zero,half,one,two,three/0.0d0,0.5d0,1.0d0,2.0d0,3.0d0/

      data bohr,cutoff/0.529177d0,1.0d+35/

      sqrt3=sqrt(three)

      rewind 66

      n=0

      loop=4
      if(iptsflag.eq.0) loop=1

      do 1001 i=1,loop

      read(66) nsubpoints

      call fastrd(66,esppt(n+1,1),nsubpoints)
      call fastrd(66,esppt(n+1,2),nsubpoints)
      call fastrd(66,esppt(n+1,3),nsubpoints)

 1001 n=n+nsubpoints

      do 1002 j=1,ncentr
      do 1002 i=1,npoints
 1002 qqd(i,j)=one/sqrt((esppt(i,1)-coord(1,j))**2
     &                 +(esppt(i,2)-coord(2,j))**2
     &                 +(esppt(i,3)-coord(3,j))**2)

      do 1003 i=1,ncontd
      do 1003 j=1,nconts

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshels(j))
      i4=i2+i3

      density(i4)=density(i4)/sqrt3
      density(i4+i1*3+3)=density(i4+i1*3+3)/sqrt3
      density(i4+i1*5+10)=density(i4+i1*5+10)/sqrt3

 1003 continue

      do 1004 i=1,ncontd
      do 1004 j=1,ncontp

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshelp(j))
      i4=i2+i3

      density(i4)=density(i4)/sqrt3
      density(i4+1)=density(i4+1)/sqrt3
      density(i4+2)=density(i4+2)/sqrt3

      density(i4+i1*3+3)=density(i4+i1*3+3)/sqrt3
      density(i4+i1*3+4)=density(i4+i1*3+4)/sqrt3
      density(i4+i1*3+5)=density(i4+i1*3+5)/sqrt3

      density(i4+i1*5+10)=density(i4+i1*5+10)/sqrt3
      density(i4+i1*5+11)=density(i4+i1*5+11)/sqrt3
      density(i4+i1*5+12)=density(i4+i1*5+12)/sqrt3

 1004 continue

      do 1005 i=1,ncontd
      do 1005 j=1,i

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nsheld(j))
      i4=i2+i3

      density(i4)=density(i4)/three
      if(i.ne.j) density(i4+1)=density(i4+1)/sqrt3
      if(i.ne.j) density(i4+2)=density(i4+2)/sqrt3
      if(i.ne.j) density(i4+3)=density(i4+3)/three
      if(i.ne.j) density(i4+4)=density(i4+4)/sqrt3
      if(i.ne.j) density(i4+5)=density(i4+5)/three

      density(i4+i1)=density(i4+i1)/sqrt3
      if(i.ne.j) density(i4+i1+3)=density(i4+i1+3)/sqrt3
      if(i.ne.j) density(i4+i1+5)=density(i4+i1+5)/sqrt3

      density(i4+i1*2+1)=density(i4+i1*2+1)/sqrt3
      if(i.ne.j) density(i4+i1*2+4)=density(i4+i1*2+4)/sqrt3
      if(i.ne.j) density(i4+i1*2+6)=density(i4+i1*2+6)/sqrt3

      density(i4+i1*3+3)=density(i4+i1*3+3)/three
      density(i4+i1*3+4)=density(i4+i1*3+4)/sqrt3
      density(i4+i1*3+5)=density(i4+i1*3+5)/sqrt3
      density(i4+i1*3+6)=density(i4+i1*3+6)/three
      if(i.ne.j) density(i4+i1*3+7)=density(i4+i1*3+7)/sqrt3
      if(i.ne.j) density(i4+i1*3+8)=density(i4+i1*3+8)/three

      density(i4+i1*4+6)=density(i4+i1*4+6)/sqrt3
      density(i4+i1*4+9)=density(i4+i1*4+9)/sqrt3
      if(i.ne.j) density(i4+i1*4+11)=density(i4+i1*4+11)/sqrt3

      density(i4+i1*5+10)=density(i4+i1*5+10)/three
      density(i4+i1*5+11)=density(i4+i1*5+11)/sqrt3
      density(i4+i1*5+12)=density(i4+i1*5+12)/sqrt3
      density(i4+i1*5+13)=density(i4+i1*5+13)/three
      density(i4+i1*5+14)=density(i4+i1*5+14)/sqrt3
      density(i4+i1*5+15)=density(i4+i1*5+15)/three

 1005 continue

      ncposs=nconts*(nconts+1)/2
      ncpopp=ncontp*(ncontp+1)/2
      ncpodd=ncontd*(ncontd+1)/2

      ncpops=ncontp*nconts
      ncpods=ncontd*nconts
      ncpodp=ncontd*ncontp

      npposs=0
      nppops=0
      nppopp=0
      nppods=0
      nppodp=0
      nppodd=0

      do 1006 i=1,nconts
      do 1006 j=1,i
 1006 npposs=npposs+nprimvc(nshels(i))*nprimvc(nshels(j))

      do 1007 i=1,ncontp
      do 1007 j=1,nconts
 1007 nppops=nppops+nprimvc(nshelp(i))*nprimvc(nshels(j))

      do 1008 i=1,ncontp
      do 1008 j=1,i
 1008 nppopp=nppopp+nprimvc(nshelp(i))*nprimvc(nshelp(j))

      do 1009 i=1,ncontd
      do 1009 j=1,nconts
 1009 nppods=nppods+nprimvc(nsheld(i))*nprimvc(nshels(j))

      do 1010 i=1,ncontd
      do 1010 j=1,ncontp
 1010 nppodp=nppodp+nprimvc(nsheld(i))*nprimvc(nshelp(j))

      do 1011 i=1,ncontd
      do 1011 j=1,i
 1011 nppodd=nppodd+nprimvc(nsheld(i))*nprimvc(nsheld(j))

      do 1012 i=1,npoints
      esppt(i,4)=zero
      do 1012 j=1,ncentr
 1012 esppt(i,4)=esppt(i,4)+charge(j)*qqd(i,j)

      do 1013 k=1,3
      do 1013 i=1,npoints
      esppt(i,k+4)=zero
      do 1013 j=1,ncentr
 1013 esppt(i,k+4)=esppt(i,k+4)
     &            -charge(j)*(esppt(i,k)-coord(k,j))*(qqd(i,j)**3)

      do 1014 i=1,npoints

      ja=1
      jb=ja+npposs
      jc=jb+ncposs

      if(jc.gt.ninteger) stop 'mespmm ss integer overflow'

      je=1
      jf=je+npposs
      jg=jf+npposs
      jh=jg+npposs
      ji=jh+npposs
      jj=ji+npposs
      jk=jj+npposs
      jl=jk+npposs
      jm=jl+npposs
      jn=jm+npposs
      jo=jn+npposs
      jp=jo+npposs
      jq=jp+npposs

      if(jq.gt.nreal) stop 'mespmm ss real overflow'

      call setupdcss(ncposs,n,nconts,nshels,mtloca,iprimvc,nprimvc,
     &               icentvc,iwkvec(ja),iwkvec(jb),cutoff,coord,alpha,
     &               coeff,wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh),
     &               wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),wkvec(jm),
     &               wkvec(jn),wkvec(jo),wkvec(jp),density)

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n

      ka=jq+1
      kb=ka+n
      kc=kb+ncposs
      kd=kc+n
      ke=kd+n
      kf=ke+n
      kg=kf+n
      kh=kg+n
      ki=kh+n
      kj=ki+n
      kk=kj+n
      kl=kk+ncposs
      km=kl+ncposs
      kn=km+ncposs

      if(kn.gt.nreal) stop 'mespmm ss real overflow'

      call efiss(ncposs,n,iwkvec(ja),iwkvec(jb),esppt(i,4),esppt(i,5),
     &           esppt(i,6),esppt(i,7),esppt(i,1),esppt(i,2),esppt(i,3),
     &           g,density,wkvec(je),wkvec(jn),wkvec(jo),wkvec(jp),
     &           wkvec(jm),wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),
     &           wkvec(ke),wkvec(kf),wkvec(kg),wkvec(kh),wkvec(ki),
     &           wkvec(kj),wkvec(kk),wkvec(kl),wkvec(km))

 1014 continue

      do 1015 i=1,npoints

      ja=1
      jb=ja+nppops
      jc=jb+ncpops*3

      if(jc.gt.ninteger) stop 'mespmm ps integer overflow'

      je=1
      jf=je+nppops
      jg=jf+nppops
      jh=jg+nppops
      ji=jh+nppops
      jj=ji+nppops
      jk=jj+nppops
      jl=jk+nppops
      jm=jl+nppops
      jn=jm+nppops
      jo=jn+nppops
      jp=jo+nppops
      jq=jp+nppops
      jr=jq+nppops
      js=jr+nppops
      jt=js+nppops

      if(jt.gt.nreal) stop 'mespmm ps real overflow'

      call setupdcps(ncpops,n,ncontp,nconts,nshelp,nshels,mtloca,
     &               iprimvc,nprimvc,icentvc,iwkvec(ja),iwkvec(jb),
     &               cutoff,coord,alpha,coeff,wkvec(je),wkvec(jf),
     &               wkvec(jg),wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),
     &               wkvec(jl),wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),
     &               wkvec(jq),wkvec(jr),wkvec(js),density)

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n

      ka=jt+1
      kb=ka+n
      kc=kb+n
      kd=kc+n*3
      ke=kd+ncpops*3
      kf=ke+n
      kg=kf+n
      kh=kg+n
      ki=kh+n
      kj=ki+n
      kk=kj+n
      kl=kk+n
      km=kl+n
      kn=km+n
      ko=kn+n
      kp=ko+n
      kq=kp+n*3
      kr=kq+n*3
      ks=kr+n*3
      kt=ks+ncpops*3
      ku=kt+ncpops*3
      kv=ku+ncpops*3

      if(kv.gt.nreal) stop 'mespmm ps real overflow'

      call efips(ncpops,n,iwkvec(ja),iwkvec(jb),esppt(i,4),esppt(i,5),
     &           esppt(i,6),esppt(i,7),esppt(i,1),esppt(i,2),esppt(i,3),
     &           g,density,wkvec(je),wkvec(jn),wkvec(jo),wkvec(jp),
     &           wkvec(jq),wkvec(jr),wkvec(js),wkvec(jm),wkvec(ka),
     &           wkvec(kb),wkvec(kc),wkvec(kd),wkvec(ke),wkvec(kf),
     &           wkvec(kg),wkvec(kh),wkvec(ki),wkvec(kj),wkvec(kk),
     &           wkvec(kl),wkvec(km),wkvec(kn),wkvec(ko),wkvec(kp),
     &           wkvec(kq),wkvec(kr),wkvec(ks),wkvec(kt),wkvec(ku))

 1015 continue

      do 1016 i=1,npoints

      ja=1
      jb=ja+ncpopp
      jc=jb+nppopp
      jd=jc+ncpopp*9

      if(jd.gt.ninteger) stop 'mespmm pp integer overflow'

      je=1
      jf=je+nppopp
      jg=jf+nppopp
      jh=jg+nppopp
      ji=jh+nppopp
      jj=ji+nppopp
      jk=jj+nppopp
      jl=jk+nppopp
      jm=jl+nppopp
      jn=jm+nppopp
      jo=jn+nppopp
      jp=jo+nppopp
      jq=jp+nppopp
      jr=jq+nppopp
      js=jr+nppopp
      jt=js+nppopp
      ju=jt+nppopp
      jv=ju+nppopp
      jw=jv+nppopp
      jx=jw+ncpopp
      jy=jx+ncpopp
      jz=jy+ncpopp

      if(jz.gt.nreal) stop 'mespmm pp real overflow'

      call setupdcpp(ncpopp,n,ncontp,nshelp,mtloca,iprimvc,nprimvc,
     &               icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),cutoff,
     &               coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &               wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &               wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &               wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &               wkvec(jw),wkvec(jx),wkvec(jy),density)

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))
      call stringcopy(n,wkvec(15*n+1),wkvec(jt))
      call stringcopy(n,wkvec(16*n+1),wkvec(ju))
      call stringcopy(n,wkvec(17*n+1),wkvec(jv))
      call stringcopy(ncpopp,wkvec(18*n+0*ncpopp+1),wkvec(jw))
      call stringcopy(ncpopp,wkvec(18*n+1*ncpopp+1),wkvec(jx))
      call stringcopy(ncpopp,wkvec(18*n+2*ncpopp+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpopp
      jy=jx+ncpopp
      jz=jy+ncpopp

      ka=jz+1
      kb=ka+n
      kc=kb+n
      kd=kc+n
      ke=kd+n*3
      kf=ke+n*3
      kg=kf+n*3
      kh=kg+n*6
      ki=kh+n*6
      kj=ki+n*9
      kk=kj+n
      kl=kk+n
      km=kl+n
      kn=km+n
      ko=kn+n
      kp=ko+n
      kq=kp+n
      kr=kq+n
      ks=kr+n
      kt=ks+n
      ku=kt+n
      kv=ku+n
      kw=kv+n
      kx=kw+n
      ky=kx+n
      kz=ky+n
      la=kz+n*3
      lb=la+n*3
      lc=lb+n*3
      ld=lc+n*3
      le=ld+n*3
      lf=le+n*3
      lg=lf+n*6
      lh=lg+n*6
      li=lh+n*6
      lj=li+n
      lk=lj+n
      ll=lk+n
      lm=ll+n*3
      ln=lm+n*3
      lo=ln+n*3
      lp=lo+n*9
      lq=lp+n*9
      lr=lq+n*9
      ls=lr+n*6
      lt=ls+n*6
      lu=lt+n*6

      if(lu.gt.nreal) stop 'mespmm pp real overflow'

      call efipp(ncpopp,n,iwkvec(ja),iwkvec(jb),iwkvec(jc),esppt(i,4),
     &           esppt(i,5),esppt(i,6),esppt(i,7),esppt(i,1),
     &           esppt(i,2),esppt(i,3),g,density,wkvec(je),
     &           wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jn),wkvec(jo),
     &           wkvec(jp),wkvec(jq),wkvec(jr),wkvec(js),wkvec(jm),
     &           wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),wkvec(ke),
     &           wkvec(kf),wkvec(kg),wkvec(kh),wkvec(ki),wkvec(kj),
     &           wkvec(kk),wkvec(kl),wkvec(km),wkvec(kn),wkvec(ko),
     &           wkvec(kp),wkvec(kq),wkvec(kr),wkvec(ks),wkvec(kt),
     &           wkvec(ku),wkvec(kv),wkvec(kw),wkvec(kx),wkvec(ky),
     &           wkvec(kz),wkvec(la),wkvec(lb),wkvec(lc),wkvec(ld),
     &           wkvec(le),wkvec(lf),wkvec(lg),wkvec(lh),wkvec(li),
     &           wkvec(lj),wkvec(lk),wkvec(ll),wkvec(lm),wkvec(ln),
     &           wkvec(lo),wkvec(lp),wkvec(lq),wkvec(lr),wkvec(ls),
     &           wkvec(lt))

 1016 continue

      do 1017 i=1,npoints

      ja=1
      jb=ja+nppods
      jc=jb+ncpods*6

      if(jc.gt.ninteger) stop 'mespmm ds integer overflow'

      je=1
      jf=je+nppods
      jg=jf+nppods
      jh=jg+nppods
      ji=jh+nppods
      jj=ji+nppods
      jk=jj+nppods
      jl=jk+nppods
      jm=jl+nppods
      jn=jm+nppods
      jo=jn+nppods
      jp=jo+nppods
      jq=jp+nppods
      jr=jq+nppods
      js=jr+nppods
      jt=js+nppods

      if(jt.gt.nreal) stop 'mespmm ds real overflow'

      call setupdcds(ncpods,n,ncontd,nconts,nsheld,nshels,mtloca,
     &               iprimvc,nprimvc,icentvc,iwkvec(ja),iwkvec(jb),
     &               cutoff,coord,alpha,coeff,wkvec(je),wkvec(jf),
     &               wkvec(jg),wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),
     &               wkvec(jl),wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),
     &               wkvec(jq),wkvec(jr),wkvec(js),density)

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n

      ka=jt+1
      kb=ka+n
      kc=kb+n
      kd=kc+n
      ke=kd+n*3
      kf=ke+n*3
      kg=kf+n*6
      kh=kg+ncpods*6
      ki=kh+n
      kj=ki+n
      kk=kj+n
      kl=kk+n
      km=kl+n
      kn=km+n
      ko=kn+n
      kp=ko+n
      kq=kp+n
      kr=kq+n
      ks=kr+n
      kt=ks+n
      ku=kt+n
      kv=ku+n
      kw=kv+n
      kx=kw+n
      ky=kx+n
      kz=ky+n
      la=kz+n
      lb=la+n*3
      lc=lb+n*3
      ld=lc+n*3
      le=ld+n*3
      lf=le+n*3
      lg=lf+n*3
      lh=lg+n*6
      li=lh+n*6
      lj=li+n*6
      lk=lj+ncpods*6
      ll=lk+ncpods*6
      lm=ll+ncpods*6

      if(lm.gt.nreal) stop 'mespmm ds real overflow'

      call efids(ncpods,n,iwkvec(ja),iwkvec(jb),esppt(i,4),esppt(i,5),
     &           esppt(i,6),esppt(i,7),esppt(i,1),esppt(i,2),esppt(i,3),
     &           g,density,wkvec(je),wkvec(jn),wkvec(jo),wkvec(jp),
     &           wkvec(jq),wkvec(jr),wkvec(js),wkvec(jm),wkvec(ka),
     &           wkvec(kb),wkvec(kc),wkvec(kd),wkvec(ke),wkvec(kf), 
     &           wkvec(kg),wkvec(kh),wkvec(ki),wkvec(kj),wkvec(kk),
     &           wkvec(kl),wkvec(km),wkvec(kn),wkvec(ko),wkvec(kp),
     &           wkvec(kq),wkvec(kr),wkvec(ks),wkvec(kt),wkvec(ku),
     &           wkvec(kv),wkvec(kw),wkvec(kx),wkvec(ky),wkvec(kz),
     &           wkvec(la),wkvec(lb),wkvec(lc),wkvec(ld),wkvec(le),
     &           wkvec(lf),wkvec(lg),wkvec(lh),wkvec(li),wkvec(lj),
     &           wkvec(lk),wkvec(ll))

 1017 continue

      do 1018 i=1,npoints

      ja=1
      jb=ja+nppodp
      jc=jb+ncpodp*18

      if(jc.gt.ninteger) stop 'mespmm dp integer overflow'

      je=1
      jf=je+nppodp
      jg=jf+nppodp
      jh=jg+nppodp
      ji=jh+nppodp
      jj=ji+nppodp
      jk=jj+nppodp
      jl=jk+nppodp
      jm=jl+nppodp
      jn=jm+nppodp
      jo=jn+nppodp
      jp=jo+nppodp
      jq=jp+nppodp
      jr=jq+nppodp
      js=jr+nppodp
      jt=js+nppodp
      ju=jt+nppodp
      jv=ju+nppodp
      jw=jv+nppodp
      jx=jw+ncpodp
      jy=jx+ncpodp
      jz=jy+ncpodp

      if(jz.gt.nreal) stop 'mespmm dp real overflow'

      call setupdcdp(ncpodp,n,ncontd,ncontp,nsheld,nshelp,mtloca,
     &               iprimvc,nprimvc,icentvc,iwkvec(ja),iwkvec(jb),
     &               cutoff,coord,alpha,coeff,wkvec(je),wkvec(jf),
     &               wkvec(jg),wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),
     &               wkvec(jl),wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),
     &               wkvec(jq),wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),
     &               wkvec(jv),wkvec(jw),wkvec(jx),wkvec(jy),density)

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))
      call stringcopy(n,wkvec(15*n+1),wkvec(jt))
      call stringcopy(n,wkvec(16*n+1),wkvec(ju))
      call stringcopy(n,wkvec(17*n+1),wkvec(jv))
      call stringcopy(ncpodp,wkvec(18*n+0*ncpodp+1),wkvec(jw))
      call stringcopy(ncpodp,wkvec(18*n+1*ncpodp+1),wkvec(jx))
      call stringcopy(ncpodp,wkvec(18*n+2*ncpodp+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpodp
      jy=jx+ncpodp
      jz=jy+ncpodp

      ka=jz+1
      kb=ka+n
      kc=kb+n
      kd=kc+n
      ke=kd+n
      kf=ke+n*3
      kg=kf+n*3
      kh=kg+n*3
      ki=kh+n*6
      kj=ki+n*6
      kk=kj+n*6
      kl=kk+n*10
      km=kl+n*10
      kn=km+n*18
      ko=kn+n
      kp=ko+n*3
      kq=kp+n
      kr=kq+n
      ks=kr+n
      kt=ks+n
      ku=kt+n
      kv=ku+n
      kw=kv+n
      kx=kw+n
      ky=kx+n
      kz=ky+n*3
      la=kz+n*3
      lb=la+n*3
      lc=lb+n
      ld=lc+n
      le=ld+n
      lf=le+n
      lg=lf+n
      lh=lg+n
      li=lh+n
      lj=li+n
      lk=lj+n
      ll=lk+n
      lm=ll+n
      ln=lm+n
      lo=ln+n*3
      lp=lo+n*3
      lq=lp+n*3
      lr=lq+n*3
      ls=lr+n*3
      lt=ls+n*3
      lu=lt+n*3
      lv=lu+n*3
      lw=lv+n*3
      lx=lw+n*6
      ly=lx+n*6
      lz=ly+n*6
      ma=lz+n*6
      mb=ma+n*6
      mc=mb+n*6
      md=mc+n*10
      me=md+n*10
      mf=me+n*10
      mg=mf+n*6
      mh=mg+n*6
      mi=mh+n*6
      mj=mi+n*18
      mk=mj+n*18
      ml=mk+n*18
      mm=ml+n*10
      mn=mm+n*10
      mo=mn+n*10

      if(mo.gt.nreal) stop 'mespmm dp real overflow'

      call efidp(ncpodp,n,iwkvec(ja),iwkvec(jb),esppt(i,4),esppt(i,5),
     &           esppt(i,6),esppt(i,7),esppt(i,1),esppt(i,2),esppt(i,3),
     &           g,density,wkvec(je),wkvec(jw),wkvec(jx),wkvec(jy),
     &           wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),
     &           wkvec(js),wkvec(jm),wkvec(ka),wkvec(kb),wkvec(kc),
     &           wkvec(kd),wkvec(ke),wkvec(kf),wkvec(kg),wkvec(kh),
     &           wkvec(ki),wkvec(kj),wkvec(kk),wkvec(kl),wkvec(km),
     &           wkvec(kn),wkvec(ko),wkvec(kp),wkvec(kq),wkvec(kr),
     &           wkvec(ks),wkvec(kt),wkvec(ku),wkvec(kv),wkvec(kw),
     &           wkvec(kx),wkvec(ky),wkvec(kz),wkvec(la),wkvec(lb),
     &           wkvec(lc),wkvec(ld),wkvec(le),wkvec(lf),wkvec(lg),
     &           wkvec(lh),wkvec(li),wkvec(lj),wkvec(lk),wkvec(ll),
     &           wkvec(lm),wkvec(ln),wkvec(lo),wkvec(lp),wkvec(lq),
     &           wkvec(lr),wkvec(ls),wkvec(lt),wkvec(lu),wkvec(lv),
     &           wkvec(lw),wkvec(lx),wkvec(ly),wkvec(lz),wkvec(ma),
     &           wkvec(mb),wkvec(mc),wkvec(md),wkvec(me),wkvec(mf),
     &           wkvec(mg),wkvec(mh),wkvec(mi),wkvec(mj),wkvec(mk),
     &           wkvec(ml),wkvec(mm),wkvec(mn))

 1018 continue

      do 1019 i=1,npoints

      ja=1
      jb=ja+ncpodd
      jc=jb+nppodd
      jd=jc+ncpodd*36

      if(jd.gt.ninteger) stop 'mespmm dd integer overflow'

      je=1
      jf=je+nppodd
      jg=jf+nppodd
      jh=jg+nppodd
      ji=jh+nppodd
      jj=ji+nppodd
      jk=jj+nppodd
      jl=jk+nppodd
      jm=jl+nppodd
      jn=jm+nppodd
      jo=jn+nppodd
      jp=jo+nppodd
      jq=jp+nppodd
      jr=jq+nppodd
      js=jr+nppodd
      jt=js+nppodd
      ju=jt+nppodd
      jv=ju+nppodd
      jw=jv+nppodd
      jx=jw+ncpodd
      jy=jx+ncpodd
      jz=jy+ncpodd

      if(jz.gt.nreal) stop 'mespmm dd real overflow'

      call setupdcdd(ncpodd,n,ncontd,nsheld,mtloca,iprimvc,nprimvc,
     &               icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),cutoff,
     &               coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &               wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &               wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &               wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &               wkvec(jw),wkvec(jx),wkvec(jy),density)

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))
      call stringcopy(n,wkvec(15*n+1),wkvec(jt))
      call stringcopy(n,wkvec(16*n+1),wkvec(ju))
      call stringcopy(n,wkvec(17*n+1),wkvec(jv))
      call stringcopy(ncpodd,wkvec(18*n+0*ncpodd+1),wkvec(jw))
      call stringcopy(ncpodd,wkvec(18*n+1*ncpodd+1),wkvec(jx))
      call stringcopy(ncpodd,wkvec(18*n+2*ncpodd+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpodd
      jy=jx+ncpodd
      jz=jy+ncpodd

      ka=jz+1
      kb=ka+n
      kc=kb+n
      kd=kc+n
      ke=kd+n
      kf=ke+n
      kg=kf+n*3
      kh=kg+n*3
      ki=kh+n*3
      kj=ki+n*3
      kk=kj+n*6
      kl=kk+n*6
      km=kl+n*6
      kn=km+n*6
      ko=kn+n*10
      kp=ko+n*10
      kq=kp+n*10
      kr=kq+n*15
      ks=kr+n*15
      kt=ks+n*18
      ku=kt+n*30
      kv=ku+n*36
      kw=kv+n
      kx=kw+n*3
      ky=kx+n*6
      kz=ky+n
      la=kz+n
      lb=la+n
      lc=lb+n
      ld=lc+n
      le=ld+n
      lf=le+n
      lg=lf+n
      lh=lg+n*3
      li=lh+n*3
      lj=li+n*3
      lk=lj+n*6
      ll=lk+n*6
      lm=ll+n*6
      ln=lm+n
      lo=ln+n
      lp=lo+n
      lq=lp+n
      lr=lq+n
      ls=lr+n
      lt=ls+n
      lu=lt+n
      lv=lu+n
      lw=lv+n
      lx=lw+n
      ly=lx+n
      lz=ly+n
      ma=lz+n
      mb=ma+n
      mc=mb+n
      md=mc+n*3
      me=md+n*3
      mf=me+n*3
      mg=mf+n*3
      mh=mg+n*3
      mi=mh+n*3
      mj=mi+n*3
      mk=mj+n*3
      ml=mk+n*3
      mm=ml+n*3
      mn=mm+n*3
      mo=mn+n*3
      mp=mo+n*6
      mq=mp+n*6
      mr=mq+n*6
      ms=mr+n*6
      mt=ms+n*6
      mu=mt+n*6
      mv=mu+n*6
      mw=mv+n*6
      mx=mw+n*6
      my=mx+n*10
      mz=my+n*10
      na=mz+n*10
      nb=na+n*10
      nc=nb+n*10
      nd=nc+n*10
      ne=nd+n*15
      nf=ne+n*15
      ng=nf+n*15
      nh=ng+n*6
      ni=nh+n*6
      nj=ni+n*6
      nk=nj+n*10
      nl=nk+n*10
      nm=nl+n*10
      nn=nm+n*15
      no=nn+n*15
      np=no+n*15
      nq=np+n*18
      nr=nq+n*18
      ns=nr+n*18
      nt=ns+n*30
      nu=nt+n*30
      nv=nu+n*30
      nw=nv+n*36
      nx=nw+n*36
      ny=nx+n*36

      if(ny.gt.nreal) stop 'mespmm dd real overflow'

      call efidd(ncpodd,n,iwkvec(ja),iwkvec(jb),iwkvec(jc),esppt(i,4),
     &           esppt(i,5),esppt(i,6),esppt(i,7),esppt(i,1),esppt(i,2),
     &           esppt(i,3),g,density,wkvec(je),wkvec(jw),wkvec(jx),
     &           wkvec(jy),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &           wkvec(jr),wkvec(js),wkvec(jm),wkvec(ka),wkvec(kb),
     &           wkvec(kc),wkvec(kd),wkvec(ke),wkvec(kf),wkvec(kg),
     &           wkvec(kh),wkvec(ki),wkvec(kj),wkvec(kk),wkvec(kl),
     &           wkvec(km),wkvec(kn),wkvec(ko),wkvec(kp),wkvec(kq),
     &           wkvec(kr),wkvec(ks),wkvec(kt),wkvec(ku),wkvec(kv),
     &           wkvec(kw),wkvec(kx),wkvec(ky),wkvec(kz),wkvec(la),
     &           wkvec(lb),wkvec(lc),wkvec(ld),wkvec(le),wkvec(lf),
     &           wkvec(lg),wkvec(lh),wkvec(li),wkvec(lj),wkvec(lk),
     &           wkvec(ll),wkvec(lm),wkvec(ln),wkvec(lo),wkvec(lp),
     &           wkvec(lq),wkvec(lr),wkvec(ls),wkvec(lt),wkvec(lu),
     &           wkvec(lv),wkvec(lw),wkvec(lx),wkvec(ly),wkvec(lz),
     &           wkvec(ma),wkvec(mb),wkvec(mc),wkvec(md),wkvec(me),
     &           wkvec(mf),wkvec(mg),wkvec(mh),wkvec(mi),wkvec(mj),
     &           wkvec(mk),wkvec(ml),wkvec(mm),wkvec(mn),wkvec(mo),
     &           wkvec(mp),wkvec(mq),wkvec(mr),wkvec(ms),wkvec(mt),
     &           wkvec(mu),wkvec(mv),wkvec(mw),wkvec(mx),wkvec(my),
     &           wkvec(mz),wkvec(na),wkvec(nb),wkvec(nc),wkvec(nd),
     &           wkvec(ne),wkvec(nf),wkvec(ng),wkvec(nh),wkvec(ni),
     &           wkvec(nj),wkvec(nk),wkvec(nl),wkvec(nm),wkvec(nn),
     &           wkvec(no),wkvec(np),wkvec(nq),wkvec(nr),wkvec(ns),
     &           wkvec(nt),wkvec(nu),wkvec(nv),wkvec(nw),wkvec(nx))

 1019 continue

      nstore=0
 
      do 1020 j=1,nunique
      do 1020 k=1,nunique
      qqa(j,k)=zero
      do 1020 i=1,npoints
 1020 qqa(j,k)=qqa(j,k)+qqd(i,k)*qqd(i,j)

      do 1021 i=1,nunique
 1021 qqb(i)=zero

      do 1022 j=1,nunique
      do 1022 i=1,npoints
 1022 qqb(j)=qqb(j)+esppt(i,4)*qqd(i,j)

      rewind 91

      do 1023 i=1,npoints
 1023 write(91,1024) esppt(i,1)*bohr,esppt(i,2)*bohr,esppt(i,3)*bohr,
     &               esppt(i,4),
     &               esppt(i,5),esppt(i,6),esppt(i,7)
 1024 format(4f20.12,/,3f20.12)

      return
      end
