integer ndim,mspin,mpart,nppss,nparts,indices,ispin,mdim,mnpair
     .,mnpss,mnpss2,mno,norbits,mngv
     .,ntypes,mnc,ncomps,nspins,nvmax,jlindem,mkmd,mnkv,mnsh,lptable
     .,mnkv2,mnshex,lgrad,lcgrad,lbuff,lbobuf,lnptbuf
     .,nkact,kmult,nvects,nvects2,ngblks,mwait,nwaits
     .,nksofk,verbose,igofr,lrsq
     .,lgen,ltf,lphase,ldet,le,lei,lv,lvee,lvep,lvpp
     .,lke,lgp2,lx,lbuffm,nwalks,ldel2u
     .,lage,ltheta,phase_dist,lweight,lwbuff,neblocks,nkmd
     .,nfpty,nlpty,idtype,nelects,ifpair,nphases,load,nstps,mngofr,lbox
     .,lvol,nprotons,idebug,inpt,mptsa,ifscl,nparms,mparms,ngvects
     .,nbact,n3act,ite,itp,ifbounce,direction,writeperiod,nwritten
     &,twist_freq,irun,nx_theta,isample_theta,fourier,nbasis

      integer mplong
      parameter (mplong=1) ! mh mplong=0: full longrange (longrange.f)
                           !    mplong>0: grad for p not evaluated (lrange.f)
      parameter (nvmax=110) ! # phases

      integer meslices,neslices,head,head_save(nvmax)  !REPTATE
      integer mpslices,npslices      !PIMC
      parameter (meslices=101,mpslices=10,mngv=400)       !REPTATE
      parameter (mparms=8) ! maximum number of searchable parameters

!mmorales
      integer nproc_world, myrank_world,maxproc_world
      integer my_comm_local,comm_world 
      parameter (maxproc_world=500)

      integer nproc,myrank,iseed,maxproc
      parameter (maxproc=110)
      integer istatmpi
      logical root,root_world

      common/mpi/nproc,myrank,iseed,root,nproc_world,myrank_world
     .,root_world,istatmpi(maxproc),my_comm_local,comm_world
     
      real*8 pi,twopi

      parameter (ndim=3,mspin=2,mpart=200,mdim=3,mwait=100,lwbuff=12
     .,mnpss=mpart/4,mnpss2=mpart/2,mno=2)
      parameter(mnkv=900,mnsh=200,lptable=400,mnkv2=mnkv*2,mnshex=900) 
      parameter (mkmd=300,mnc=2,lbobuf=11,lnptbuf=11)
      parameter (ltf=1,lage=2,ldet=3,lphase=4,le=5,lv=6,lvee=7,lvep=8
     . ,lvpp=9,lke=10,lgp2=11,ldel2u=12,lei=13,ltheta=14
     . ,lrsq=ltheta+ndim+1,lgen=lrsq+mnc+1
     . ,lweight=lgen+1,lx=lweight+lwbuff+lbobuf+lnptbuf
     . ,lbox=lweight+lwbuff+1,lvol=lbox+9
     . ,lbuffm=lx+4*ndim*mpart ,mnpair=mpart*(mpart-1)/2
     . ,mngofr=mnc+2,mptsa=500)
      parameter (pi=3.1415 92653 58979d0)
      parameter (twopi=2.d0*3.1415 92653 58979d0)

      logical debug,npt,singletwist
      character*8 pname,sparms
      character molpot*6
      common/cint/pname(mnc),sparms(8),molpot

      common/sint/nppss(mspin,mnc),nparts,indices(mpart),ispin(mpart)
     .,nkact(2),nksofk,kmult(0:mnsh,2),nvects(2)
     .,nvects2(2),ngblks,nwaits(0:mwait),nparms
     .,verbose,lgrad,lcgrad,lbuff,nkmd
     .,igofr(mnpair,2),phase_dist(0:21),neblocks(0:lwbuff+1),ncomps(mnc)
     .,ntypes,nspins(mnc),nfpty(mnc),nlpty(mnc),idtype(mpart),nelects
     .,nprotons,ifpair,nphases,load,nstps(nvmax),debug,jlindem,npt
     .,idebug,inpt,neslices,ifscl,npslices,ngvects,singletwist
     .,nbact,n3act,ite,itp,ifbounce,direction,writeperiod,nwritten
     &,norbits,twist_freq,irun,nx_theta,isample_theta
     .,fourier(mdim,nvmax),nbasis

      real*8 charge,rho,beta,hbs2m,theta,cutr,cutr2,cutr2min,enorm
     .,xold,rs,rkcomp,ulr,usr,tauc,wsites,rknorm,dtable,dtablei
     .,rhoke,rhokp,pwmate,kmod,kvec,etrial,urpa,enum,edenom,ebav,ebsq
     .,ensum,edsum,tasnorm,taseni,taseni2,tase,tasecor,tase2,wtbl
     .,phweight,totload,cut1,cut2,cut3,psites,pext,hbox,hboxi,volume
     .,vlog,dim_cutoff,pwmatp,pwmat,psit,gamma,oep,theta_new,pathwt
     .,cuta2,r1a,atablei,ept,temperaturea,ehbs2m,minoverlap
     .,ubackflow,u3back,u3jas,bondvector,bdlmx,bdlmn ,bparms,dbparms
     .,deltatwist


      common/sreal/xold(lbuffm),theta(mdim,nvmax),rkcomp(ndim,mnkv,2)
     .,ulr(0:mnshex,mnc,mnc,5,2),usr(0:lptable,4,mnc,mnc,8,2)
     .,rknorm(0:mnshex,2),dtable(2),dtablei(2),wsites(ndim,mpart)
     .,psites(mdim,mpart),rhoke(mnkv2,2),rhokp(mnkv2,mpslices,2)
     .,pwmate(mnkv2,mpart),pwmatp(mnkv2,mpart,mpslices,2)
     .,pwmat(mnkv2,mpart),psit(2,meslices),gamma,oep(4),pathwt(3)
     .,etrial,tauc,kvec(ndim,mkmd,2),kmod(mkmd,2)
     .,urpa(0:lptable,mnc,mnc,2),wtbl(0:lptable,4,mnc,2)
     .,enum(0:lwbuff+1),edenom(0:lwbuff+1),theta_new(3)
     .,edsum(0:lwbuff+1),ensum(0:lwbuff+1),ebav(0:lwbuff+1)
     .,ebsq(0:lwbuff+1),hbs2m(mnc),charge(mnc),rho,cutr(2),cutr2(2),beta
     .,enorm,cutr2min,bparms(mparms,2),dbparms(mparms)
     .,tasnorm,taseni,taseni2,tase,tasecor,tase2,phweight(nvmax),rs
     .,totload,cut1,cut2,cut3,pext,hbox(mdim,mdim,2),hboxi(mdim,mdim,2)
     .,volume(2),vlog(2),dim_cutoff,minoverlap
     .,ubackflow(0:2*lptable,mnc,mnc,2),u3back(0:2*lptable,mnc,mnc,2) !mh
     .,u3jas(0:2*lptable,mnc,mnc,2),bondvector(ndim,mpart),bdlmx,bdlmn
     .,deltatwist

!bomc
      integer maxbuff,nptable,nbgofr
      real*8 hist_ct,grsum,dpi,electron,proton,grcu,grcu2,grpsum
      real*8 pforce,vcl,pcl,action,kinen,vsc,psc,scforce,deloc,dact
      real*8 pfkin,hist_ovlp,epsilon,eref,etref,stdev,ocfcu,ocfcu2
      real*8 molcenter,bond,bond_hist,orient_hist,bhcu,ohcu,ocfsum
      real*8 grcmsum,grcmcu,grcmcu2,bhcu2,ohcu2
      parameter(maxbuff=2201,nptable=200)
      common/cbomc/proton(mdim,mpart,mpslices,2)
     .,electron(mdim,mpart,meslices,nvmax) ! REPTATE
     .,molcenter(mdim,mpart/4,mpslices,2) ! molecules
     .,bond(mdim,mpart/4,mpslices,2) ! molecules
     .,bond_hist(nptable+2),orient_hist(nptable+2)
     .,bhcu(nptable),ohcu(nptable),bhcu2(nptable),ohcu2(nptable)
     &,pforce(mdim,mpart,mpslices,2),pfkin(mdim,mpart,mpslices,2)
     &,hist_ct(maxbuff+1),hist_ovlp(nptable+1),scforce(mdim,mpart,2)
     &,grsum(nptable+2,mngofr),grcu(nptable,mngofr),pcl(2)
     &,grcu2(nptable,mngofr),dpi,vcl(2),grpsum(nptable+2)
     .,ept(0:mptsa,4,3),action(2),kinen(2),vsc(2),psc(2),deloc(2)
     .,dact(2),cuta2,r1a,atablei,temperaturea,ehbs2m,epsilon,eref(2)
     .,etref(2,nvmax),stdev(nvmax),ocfsum(nptable+2),ocfcu(nptable)
     .,ocfcu2(nptable),grcmsum(nptable+2),grcmcu(nptable)
     .,grcmcu2(nptable)
      common/cbomci/head,head_save,nbgofr

! correlation fuctions s(k),g(r)
      real*8 gofr,gsum,gsq,sofk,sksum,sksq,skcnt
      common /cfcm/skcnt,gsum(0:lptable),gsq(0:lptable)
     .  ,sksum(mnsh),sksq(mnsh),sofk(mnsh),gofr(0:lptable)

! variational parameters
       integer ifbackflow, if3body,hedf2,pot_expon,zerotwist,iself
     .,ifbacka,ifbackn,ifselfa,ifextraa,if3bodya,otype   !mh
     .,if3bodyn     !mh
     .,ifbacknee,if3bodynee  !cp
       real*8 plambdab,psb,prb,pwb,plambdat,prt,pwt,pextra,pwe,pc
     .,plambdabee,psbee,prbee,pwbee,plambdatee,prtee,pwtee,pb
     .,pce,pde,pee,pself,pws
     .,pla,pls,ple,p3ep,p3ee,pbee,pbep   !mh
     .,pa1,pa2,pb1,ps1,ps2,ps3,ps4  ! cp

       common/parameteri/ifbackflow,if3body,hedf2,pot_expon,zerotwist
     .,iself,ifbacka,ifbackn,ifselfa,ifextraa,ifbacknee,if3bodynee  
     .,if3bodya,if3bodyn,otype(mno)
       common/parameterr/pc(mno),pb(mno)
     .,plambdab,psb,prb,pwb,plambdat,prt,pwt
     .,plambdabee,psbee,prbee,pwbee,plambdatee,prtee,pwtee   
     .,pextra,pwe,pce,pde,pee,pself,pws,pla,pls,ple
     .,p3ep,p3ee,pbee,pbep !mh
     .,pa1,pa2,pb1,ps1,ps2,ps3,ps4  ! cp

! pimc
      integer mtab,mlevels,mterms,mpfile,ntermact,ntermeng,ipfilep,ifew
      integer iphase,islice,nsl,nmolecules,pinmolecules,companion,tabmem
      integer molmoves
      parameter (mtab=505,mlevels=1,mterms=11,mpfile=2)
      real*8 vtable,dmtable,uka,rpfilep,cutrp,taup,hbar2m,vel,act
     .,grad2pimc,vpimc,virpimc,kspring,kinact,taup1,ukep,grad2mix
     .,vimage,vtail
      common/rpimc/vtable(mtab,4),dmtable(mtab,4,mterms,mlevels,2)
     .  ,uka(mnkv2,mnc,mnc,3,mlevels),vel(mdim,mpart,mpslices,2)
     .  ,ukep(0:mnshex,3,2),rpfilep(6,mpfile),act(2),vpimc(2)
     .  ,virpimc(2),grad2pimc(2),kspring(2),kinact(2),grad2mix(2)
     .,cutrp,taup,hbar2m,taup1,vimage,vtail
      common/ipimc/ipfilep(4,mpfile),ntermact,ntermeng,ifew,tabmem
     .  ,iphase(maxproc),islice(mpslices,maxproc)
     .  ,nsl(maxproc),nmolecules,pinmolecules(mnpss,2),companion(mpart)
     .  ,molmoves

        real*8 drij,rij,rij_i,rij_rem,bond_length
        integer   rij_ind
        common/rdistance/drij(ndim,mpart,mpart),rij(mpart,mpart)
     & ,rij_rem(mpart,mpart),rij_i(mpart,mpart),bond_length
        common/idistance/rij_ind(mpart,mpart)


      real*8 amat,bmat,qp,orb,del2orb,gradorb,fulld2orb
      common/cbackflow/amat(ndim,ndim,mpart,mpart),qp(ndim,mpart)
     .,bmat(ndim,mpart),orb(mnpss,mnpss2),del2orb(mnpss,mnpss2)
     .,gradorb(mnpss,mnpss2,ndim),fulld2orb(mnpss,mnpss2,ndim,ndim)

! kdelaney: findbands and cslater variables
      integer mnpw,nocc,oeptype,ngvects_oep,maxphperproc,mngd,maxg
      parameter (maxphperproc = 10) !Maximum number of phases per CPU...for cmat
      parameter(mnpw = 4500) !Max number of planewaves for coefficients
      parameter(mngd = 40)   !Max number of planewaves in 1 spatial dimension
      logical useoep,usedft,oepnoupdate
      integer printbanddata !Extra screen output? 0=none, 1=basic, 2=detailed
      logical etotcomp !When doing DFT, converge w.r.t. Etot or Mu?
      logical cusp_corr !Correct e-p cusp for finite-basis?
      !Method for SCF cycle in DFT. 1=simple mix, 2=Broyden
      integer scfmethod,smeartype,pwlist
      real*8 ecut,dftmixing,mixmaxg,smearwidth,pwcart
      real*8 bt_mwl !Max. weight loss in basis truncation
      complex*16 cmat !C(G) coefficients for eigenfunctions

      !Common blocks
      common/ibands/nocc,oeptype,ngvects_oep(2,maxphperproc,mpslices),
     .              printbanddata,scfmethod,smeartype,
     .              pwlist(mdim,mnpw,2,maxphperproc,mpslices),
     .              maxg(mdim,2,maxphperproc,mpslices)

      common/rbands/ecut,dftmixing,bt_mwl,mixmaxg,smearwidth,
     .              pwcart(mdim,mnpw,2,maxphperproc,mpslices)
      common/cbands/cmat(mnpss,mnpw,2,maxphperproc,mpslices) !2 for different configs
      common/lbands/useoep,usedft,etotcomp,cusp_corr,oepnoupdate

      !Which solver method will be used?
      !applyhwithfft: true => use FFT method for computing h|evec> in iterdiag
      !             : false => use BLAS multiply
      !useriterdiag : true => Use iterative diagonaliser for lowest nocc eigenstates
      !             : false => use LAPACK solver for lowest nocc eigenstates
      logical applyhwithfft,useiterdiag
      common/lbsolver/applyhwithfft,useiterdiag

!     Logical variables for enable/disable individual RPA jastrows
      logical eejastrow,epjastrow
      common/ljastrow/eejastrow,epjastrow
      
!     Variables for keeping the p-p, p-e and e-e potential energy
!     Separate. These are used to add up over particles in pairtf and
!     lrange. They are NOT everages over electronic configurations.
!     They are initialised on each cycle and put into xnew (in calpps) to
!     be averaged in reptate and impsamp2. They are then used to pass the data to add_av

      real*8 peee,pepp,peep
      real*8 ppeee,ppeep,ppepp !inow-dependent for impsamp averaging
      common/pecompscmn/peee,pepp,peep,ppeee(2),ppeep(2),ppepp(2) 
! kdelaney: end changes

! mmorales: begin changes
! to study relaxation properties of RQMC during EDIFF calculation
      integer ediff_oep,save_warmup,pstep_warmup,pstep_w,first_getc
      integer nblq_warmup,getconf_warmup,num_warmup_blq,getc_type 
      integer oneside_pt,oneside  

      common/ediffFix/ediff_oep,save_warmup,pstep_warmup,pstep_w
     .,nblq_warmup,getconf_warmup,num_warmup_blq,getc_type,first_getc
     .,oneside_pt,oneside 
! mmorales: end changes