! Born-Oppheneimer MC for Hydrogen
! bopimc3.09: BOMC with NPT-REPTATION-PIMC-TABC 13/3/2003
!             fully numerical and fully analitycal twfs
!	      old and new way of selecting phases
! 
      subroutine mctas(qid_local) 
      implicit none
      integer mnprm,mfreq
      parameter (mnprm=20,mfreq=4)

      include 'tas.cm'
      include 'caver.cm'
      include 'jpoint.cm'

      character qid*20,fs*20,qtime*26,p(mnprm)*16,qid_local*20
     .,name*10,frestart*15,punits(14)*10,padd(mnprm)*16
      logical ife
      integer nexthit(mfreq),restart,ncell(ndim),freq(mfreq)
     .,ifwarm,iwarm,la,intread,ln,ie,ip,ind,nadd,ikind_th,je
     .,nxtalp,it,ifind,k,n,ipickoff,lpx,lpx1,l,i,j,is,ik,nptsa,jp
      real*8 rnd,rhoi,rlread,timer,time_initial,pbb(3),dum
     .,rsc(ndim,mpart),stemp(mdim,mpart)
     .,rsa,rna,dummy,cuta,dparms,xxx
cpierleo
      integer ithet
      real*8 rng,rng_cm
      external rng,rng_cm
cpierleo
      include 'cpbc.cm'


cpierleo MPI
!      call parallel_init()
      call second(time_initial)
      if (nproc.gt.maxproc) then
       if (root) write (6,*) 'nproc.gt.maxproc',maxproc
       call parallel_end()
      endif
!      if (root) then
!       write (6,*) ' input run id (lt 12 characters)'
!      write (6,*) ' using bopimc.in by default.'
!       read (*,'(a)') qid
      qid = qid_local
!      endif
! broadcast qid among processors
!      call bcast_char(qid,20,0)
      lpx=index(qid,' ')-1
!     write (6,'(i5,a20,i3)') myrank,qid,lpx
      if (root) then
       inquire(file=qid(1:lpx)//'.out',exist=ife)
       open(6,file=qid(1:lpx)//'.out',form='formatted',status='unknown')
       if(ife) then
         write (6,*)' output file already exists'
       endif
       write (6,*) 
     &  'BOMC for hydrogen with NPT+REPT+TABC+PIMC 8/23/02'
       write (6,*) ' run id = ',qid
       call timedate(root,qtime)
!kdelaney: open files for CASTEP input files. Saves INITIAL geometry
!       fs = qid(1:lpx)//'.castep.cell'
!       open(35,file=fs,status='unknown')
!       fs = qid(1:lpx)//'.castep.param'
!       open(36,file=fs,status='unknown')
!kdelaney: end modifications
      endif
cpierleo

      fs=qid(1:lpx)//'.in'
      if (root) write (6,*) ' trying to open file ',fs
      open(1,file=fs,status='old',form='formatted')


! DEFINE DEFAULTS HERE
!     temperature=0.d0
      restart=0
      neslices=1
      npslices=1
      norbits=0
      head=1
      nkmd=0
      ikind_th=1
      isample_theta=0
      ifpair=1
      ifbackflow=0
      if3body=0
      hedf2=0
      pot_expon=0
      cuta2=0.d0 ! default is no quantum action
      cuta=0.d0
      writeperiod=0
      nparms=0
      twist_freq=0
      nmolecules=0
      molmoves=0
      irun=2

! kdelaney: band code defaults
      ecut=0.d0
      useoep = .false.
      oepnoupdate = .false.
      usedft = .false.
      dftmixing = 0.5d0
!     Which OEP diagonalizer method to use? Default = LAPACK
      useiterdiag = .false.
      applyhwithfft = .false.
!     Jastrow defaults (all on)
      eejastrow = .true.
      epjastrow = .true.
      cusp_corr = .true.
! kdelaney: end changes

! default for cutoff of ewald sums  (cutr*cutk=dim_cutoff)
      if(ndim.eq.2) then 
         dim_cutoff=23.56
      else
         dim_cutoff=18.849556d0
      endif

      verbose=0
      pextra=0.d0
      pwe=0.d0
      pself=0.d0
      pws=0.d0
      do l=1,mno
       pc(l)=0.d0
       pb(l)=0.d0
      enddo
      cutr(1)=0.d0
      do n=1,nvmax
       phweight(n)=0.d0
       do l=1,ndim
        theta(l,n)=0.d0
       enddo
      enddo
      do l=1,ndim
       ncell(l)=0
      enddo
      zerotwist=1
      nkact(1)=0
      nksofk=0
      timer=1.d100
      minoverlap=0.d0
      ntypes=0
      nparts=0
      iseed=777
      jlindem=0
      debug=.false. !.true.
      idebug=0
      do la=1,lbuffm
       xold(la)=0.d0
      enddo
! cp initialization if the variational parameters
      plambdab=0.d0
      prb=0.d0
      pwb=0.d0
      plambdat=0.d0
      prt=0.d0
      pwt=0.d0
      plambdabee=0.d0
      prbee=0.d0
      pwbee=0.d0
      plambdatee=0.d0
      prtee=0.d0
      pwtee=0.d0
      pa1=0.d0
      pa2=0.d0
      pb1=1.d0
      ps1=0.d0
      ps2=0.d0
      ps3=0.0d0
      ps4=0.d0
      tabmem=1
! cp end      
      epsilon=1.0d0
      eref(1)=0.d0
      eref(2)=0.d0
! mh initialization for analytical functions
      ifbacka=0
      ifbackn=0
      ifbacknee=0
      pla=1.d0
      ifselfa=0
      pls=1.d0
      ifextraa=0
      ple=1.d0
      if3bodya=0
      if3bodyn=0
      if3bodynee=0
      p3ee=0.d0
      p3ep=1.d0
      nbact=0
      n3act=0
! end mh

!mmorales
      save_warmup=0
      ediff_oep=0
      pstep_warmup=0 
      pstep_w=1 
      nblq_warmup=2 
      getconf_warmup=0
      num_warmup_blq=2
      getc_type=1
      first_getc=0

! random number initialization with seed extraction
      call easy_init_rng(iseed,myrank_world,nproc_world,idebug)
      debug=.false.
      idebug=0
      npt=.false.
      singletwist=.false.

!read input 2 times: for system info and for running.
1     j=ipickoff(1,p,n,mnprm,root)
       if(j.ne.0) go to 11
       if(n.le.0) go to 1

       if(p(1).eq.'PARTICLE') then  ! add a new type of particle
            if(ntypes.eq.mnc) then
              if (root) write (6,*)' mnc too small ',ntypes+1,mnc
              call parallel_end()
            endif
            ntypes=ntypes+1
            pname(ntypes)=p(2)

            nfpty(ntypes)=nparts+1
            ncomps(ntypes)=intread(p(3))
            nparts=nparts+ncomps(ntypes)
            nlpty(ntypes)=nparts
            do i=nfpty(ntypes),nlpty(ntypes)
              idtype(i)=ntypes
            enddo

            if(ncomps(ntypes).le.0) then
             if (root) write (6,*) 'ncomps(',ntypes,')<0'
             call parallel_end()
            endif
            if(nparts.gt.mpart) then
             if (root) write (6,*) 'nparts.gt.mpart'
             call parallel_end()
            endif
 
            hbs2m(ntypes)=rlread(p(4))
            if(hbs2m(ntypes).lt.0.d0)then 
             if (root) write (6,*) 'hbs2m(',ntypes,')<0'
             call parallel_end()
            endif
            charge(ntypes)=rlread(p(5))
            nspins(ntypes)= intread(p(6))
            if(nspins(ntypes).gt.mspin) then
             if (root) write (6,*) 'nspins(ntypes).gt.mspin'
             call parallel_end()
            endif
    
            if(nspins(ntypes).gt.0) then
             nppss(nspins(ntypes),ntypes)=ncomps(ntypes)
             do k=1,nspins(ntypes)-1
              nppss(k,ntypes)=intread(p(6+k))
              nppss(nspins(ntypes),ntypes)=nppss(nspins(ntypes),ntypes)
     &                                    -nppss(k,ntypes)
             enddo
             if (root) write (6,*)' spins ',nspins(ntypes)
     &                 ,(nppss(k,ntypes),k=1,nspins(ntypes))
             do k=1,nspins(ntypes) ! number of single particle states needed
              nkmd=max(nkmd,nppss(k,ntypes))
             enddo
            endif

       elseif(p(1).eq.'LATTICE') then ! set initial conditions of the particles
            it=ifind(p(2),pname,ntypes)
            nxtalp=0  
            rhoi=ncomps(it)/volume(1)
           if(n.ge.4) then !extra parameters for molecular lattice
               nxtalp=intread(p(4))
               if(n.ge.4+ndim) then
                 do l=1,ndim
                  ncell(l)=intread(p(4+l))
                 enddo
                 do l=5+ndim,n
                 wsites(l-4-ndim,1)=rlread(p(l))
                 enddo
               else
                 wsites(1,1)=0.d0
               endif
           else
               nxtalp=0
                do l=1,ndim
                 ncell(l)=0
                enddo
                wsites(1,1)=0.d0
           endif
           call sites(wsites,ncomps(it),ndim,nxtalp,rhoi,ncell,ell,root) 
           do l=1,ndim
              elli(l) = 1.0d0/ell(l)
           enddo
            if(n.gt.2) then
                rnd=2*rlread(p(3))
            else
                rnd=0.d0
            endif

            if(p(2).eq.'p') then
! find molecules : it is necessary THAT THE CARD "MOLECULES" HAS BEEN ALREADY ISSUED
             if(nmolecules.eq.0) then
              ind=nfpty(it)-1
              do i=1,ncomps(it)
               do l=1,ndim
                rsc(l,ind+i)=wsites(l,i)+rnd*(rng_cm()-.5)
                psites(l,i)=wsites(l,i)
               enddo
              enddo
             elseif(nmolecules.lt.0) then
              call setup_molecules
              do k=1,nmolecules
               i=pinmolecules(k,1)-nfpty(it)+1
               j=pinmolecules(k,2)-nfpty(it)+1
               do l=1,ndim
                xxx=rnd*(rng_cm()-.5)
                rsc(l,i)=wsites(l,i)+xxx
                rsc(l,j)=wsites(l,j)+xxx
                psites(l,k)=molcenter(l,k,1,1)
               enddo
              enddo
             else
              write (*,*) 'mctas: nmolecules incorrect:',nmolecules
              call parallel_end()
             endif
            else
             ind=nfpty(it)-1
             do i=1,ncomps(it)
              do l=1,ndim
               rsc(l,ind+i)=wsites(l,i)+rnd*(rng()-.5)
              enddo
             enddo
            endif

       elseif(p(1).eq.'NESLICES') then
            neslices=intread(p(2))
            if(neslices.le.0.or.neslices.gt.meslices) then
             if(root) write (6,*) 'neslices out of range '
     .        ,neslices,meslices
             call parallel_end()
            endif
            gamma=0.d0
            ifbounce=1
            if(n.gt.2) then
                gamma=rlread(p(3))
                if(gamma.lt.0.d0.or.gamma.gt.1.d0-2.d0/neslices) then
                 if(root)write (6,*)'gamma out of range ',gamma,neslices
                 call parallel_end()
                endif
                if(n.gt.3) then
                ifbounce=intread(p(4))
                endif
            endif
           direction=1
           if(root)write(6,*) 'number of electron slices=',neslices
     .     ,' interior prob=',gamma,' bounce=',ifbounce
       elseif(p(1).eq.'SEMICL') then
            name=p(3)
            lpx1=index(name,' ')-1
!            name=name(1:lpx1)//'.dmm'
            inquire(file=name(1:lpx1)//'.dmm',exist=ife)
            if(ife) then
             open(3,file=name(1:lpx1)//'.dmm',form='formatted')
            else
             if(root) write(6,*) 'file dmm does not exists'
             call parallel_end()
            endif
            read(3,*) ehbs2m,rsa,temperaturea,r1a,rna,nptsa
            if(nptsa.gt.mptsa)then
             if(root) write(6,*) 'mctas: nptsa.gt.mptsa'
             call parallel_end()
            endif
            it=ifind(p(2),pname,ntypes)
!           if (ehbs2m.ne.2.d0*hbs2m(ip)) then
!            if(root) then
!             write(6,*)'ehbs2m.ne.2.d0*hbs2m(ip)'
!             write(6,*)'ehbs2m=',ehbs2m,' 2.d0*hbs2m(ip)=',2.*hbs2m(ip)
!            endif
!            call parallel_end()
!           endif
            if(n.ge.4) then
             cuta=rlread(p(4))
             cuta2=cuta**2
            else
             cuta2=rna**2
            endif
            atablei=(nptsa-1.d0)/(rna-r1a)
            do i=0,nptsa-1
             read (3,*)dummy,ept(i,1,1),ept(i,1,2),ept(i,1,3)
            enddo
            do l=1,3
             call spline(ept(0,1,l),nptsa,mptsa,1)
            enddo
            close(3)
            if(n.ge.5) ifscl=intread(p(5))
       elseif(p(1).eq.'NPSLICES') then
            npslices=intread(p(2))
            if (npslices.gt.mpslices) then
             if (root) write(6,*) 'npslices.gt.mpslices'
             call parallel_end()
            endif
            nadd=n
            do i=1,n
             padd(i)=p(i)
            enddo
       elseif(p(1).eq.'WRITE_SOFK'.and.n.gt.2)then
            nksofk=max(nksofk,intread(p(3)))
            if (root) write (6,*)'number shells for s(k)=',nksofk
      elseif(p(1).eq.'WRITE_SEARCH')then
             writeperiod=intread(p(2))
             nwritten=0
             if(writeperiod.gt.0) then
             write (*,*)' writeperiod ',writeperiod
             fs=qid(1:lpx)//'.srh'
             open(45,file=fs,form='unformatted',status='unknown')
             endif
       elseif(p(1).eq.'NOPAIR') then
            ifpair=0
       elseif(p(1).eq.'ETRIAL') then
            etrial=rlread(p(2))
!      elseif(p(1).eq.'TEMPERATURE') then
!           temperature=rlread(p(2))
!           if(temp.lt.0.d0)then
!            if (root) write (6,*) 'temp.le.0'
!            call parallel_end()
!           endif
       elseif(p(1).eq.'SEED') then
            iseed=intread(p(2))
            if(root)write(6,*) 'initialize random numbers'
            call easy_init_rng(iseed,myrank_world,nproc_world,idebug)
            if (debug) write(70+myrank,*)'first Random#:',rng()
            if (debug) write(70+myrank,*)'first CM_Random#:',rng_cm()
       elseif(p(1).eq.'RESTART') then
            restart=2
            frestart=qid(1:lpx)//'.chk'
       elseif(p(1).eq.'LINDEM') then
            jlindem=1
            it=ifind('p',pname,ntypes)
            if (n.ge.2) then
             ln=index(p(2),' ')-1
             fs=p(2)(1:ln)//'.sites'
            else
             fs=qid(1:lpx)//'.sites'
            endif
            inquire(file=fs,exist=ife)
            if(ife) then
             open(8,file=fs,form='formatted',status='old')
             if (nmolecules.eq.0) then
              do i=1,ncomps(it)
               read(8,'(3g20.10)') (psites(l,i),l=1,ndim)
              enddo
             else
              do i=1,nmolecules
               read(8,'(3g20.10)') (psites(l,i),l=1,ndim)
              enddo
             endif
             if(root) write(6,*)'Lattice Sites read from file ',fs
            else
             if (root) then
              open(8,file=fs,form='formatted',status='new')
              if (nmolecules.eq.0) then
               do i=1,ncomps(it)
                write(8,'(3g20.10)') (psites(l,i),l=1,ndim)
               enddo
              else
               do i=1,nmolecules
                write(8,'(3g20.10)') (psites(l,i),l=1,ndim)
               enddo
              endif
              write(6,*)'Lattice Sites written to file ',fs
             endif
            endif
            close(8)
!           if (debug) then
!            do i=1,ncomps(it)
!             if(root)write (6,'(i5,3f10.5)') i,(psites(l,i),l=1,ndim)
!            enddo
!           endif
      elseif(p(1).eq.'TIMER') then
            timer=rlread(p(2))
       elseif(p(1).eq.'RS') then 
              rs=rlread(p(2))
              if(n.ge.3)cutr(1)=rlread(p(3))
!             if(n.ge.3)nkact(1)=intread(p(3))
!             if(nksofk.eq.0) nksofk=nkact(1) ! record s(k) to same values
              ie=ifind('e',pname,ntypes)
              ip=ifind('p',pname,ntypes)
              ite=ifind('e',pname,ntypes)
              itp=ifind('p',pname,ntypes)
              if (root) write (6,*) ' electron type is ',ie
              if (root) write (6,*) ' proton type is ',ip
              rho=ndim/(2.0d0*(ndim-1)*pi*rs**ndim)
              nelects=ncomps(ie)
              nprotons=ncomps(ip)
              enorm=nelects+nprotons
!             enorm=nelects !DMC  normalization 
              volume(1)=nelects/rho
              volume(1)=nprotons/rho
              do l=1,ndim
                ell(l)=volume(1)**(1.d0/ndim)
              enddo
       elseif(p(1).eq.'DIM_CUTOFF') then
              dim_cutoff=rlread(p(2))
       elseif(p(1).eq.'READ_STATE') then
          restart=1
          if(n.gt.1) then
             ln=index(p(2),' ')-1
             frestart=p(2)(1:ln)//'.chk'
          else
             frestart=qid(1:lpx)//'.chk'
          endif
       elseif(p(1).eq.'THETA') then
          if(intread(p(3)).ne.ndim) then
           if(root)write(6,*) 'Dimension does not match with THETA'
           call parallel_end()
          endif
          if(root)open(95,file=qid(1:lpx)//'.phload',form='formatted'
     &        ,status='unknown')
          if(p(2).eq.'SINGLE') then !single phase run
             singletwist=.true.
             nphases=1
             totload=1.d0
             do i=1,ndim
                theta(i,1)=rlread(p(i+3))
             enddo
             phweight(1)=1.d0/dble(nproc) 
             if(root) then
                do i=1,nproc
                   iphase(i)=1
                enddo
             endif
          else                  ! parallel run
             nx_theta=intread(p(3+1)) ! number of angles in x direction
                                ! and half in the y,z directions
             if(p(2).eq.'INVERSION') then ! with inversion symmetry only
               ikind_th=0
               if(n.ge.5)then
                 isample_theta=intread(p(3+2)) ! sampling if it is 1
                 if (isample_theta==1) then
                   deltatwist=rlread(p(6))
                   if(deltatwist.le.0.d0.or. deltatwist.ge.1.d0) then
                     if(root)write(*,*)
     .                 'mctas: deltatwist not appropriate:',deltatwist
                     call parallel_end()
                   endif
                   if(root) write(*,*) 
     .                 'SAMPLING TWISTS: deltatwist=',deltatwist
               endif
               endif
             elseif(p(2).eq.'EG') then
                ikind_th=1
             else
                if(root) write(*,*) 'unknown symmetry for twist phases'
                call parallel_end()
             endif
!     construct the phase table and assign a phase and a weight to any processor
             call thetaset(nx_theta,ikind_th)
          endif
       elseif(p(1).eq.'THETANEW') then
         twist_freq=intread(p(2))
         do i=1,ndim
          theta_new(i)=rlread(p(i+2))
         enddo
       elseif(p(1).eq.'EXTRA') then
           pextra=rlread(p(2))
           if(n.ge.3)pwe=rlread(p(3))
           if(n.ge.4)pce=rlread(p(4))
           if(n.ge.5)pde=rlread(p(5))
           if(n.ge.6)pee=rlread(p(6))
           if(root)write (6,*)' using extra term for e-p jastrow'
       elseif(p(1).eq.'SELF') then
           pself=rlread(p(2))
           if(n.ge.3)pws=rlread(p(3))
           iself=2*ifind('e',pname,ntypes)
           if(root)write (6,*)' using extra term for e-e jastrow',iself
       elseif(p(1).eq.'CLPOT') then
           if (p(2).eq.'YUKAWA') then
! classical prerejection potential: single or double yukawa
            pa1=rlread(p(3))  ! 1st screening parameter
            if(n.ge.4) then
             pa2=rlread(p(4)) ! 2nd screening parameter
             pb1=rlread(p(5)) ! relative weight
            endif
           else
            if(root) write(*,*) 'mctas: wrong setting for CLASSICAL'
            call parallel_end()
           endif
       elseif(p(1).eq.'CLMOLPOT') then
           if(root) write(*,*) '#molecules',nmolecules
           if (nmolecules.eq.0) then
            if(root) write(*,*)'wrong setting for CLMOLPOT'
            call parallel_end()
           endif
           molpot=p(2)
           if (p(2).eq.'YUKAWA') then
            ps1=rlread(p(3)) ! intramolecular screening
            ps2=rlread(p(4)) ! spring constant 
            ps3=rlread(p(5)) ! equilibrium distance
            ps4=rlread(p(6)) ! constant shift
           elseif (p(2).eq.'SG') then
           elseif (p(2).eq.'KW') then
           else
            write(*,*) 'mctas: wrong setting for CLMOLPOT'
            call parallel_end()
           endif
       elseif(p(1).eq.'SSEARCH') then
         nparms=n-2
         if(nparms.gt.mparms) then
          write(*,*) 'mctas:nparms.gt.mparms'
          call parallel_end()
         endif
         dparms=rlread(p(2))
         if(root)write (6,*)'number of parameters ',nparms
         do i=1,nparms
           sparms(i)=p(i+2)
           call setparm(sparms(i),bparms(i,1),1)
           dbparms(i)=dparms*abs(bparms(i,1))
         enddo
           open(41,file=qid(1:lpx)//'.par'
     &      ,form='formatted',status='unknown')
! kdelaney: New options for band-based trial functions
       elseif(p(1).eq.'OEP') then ! band potential
          if(ndim.ne.3 .and. root)then
             write(6,*)'ERROR: Findbands currently only for 3D'
             call parallel_end()
          endif
          ecut = rlread(p(2))
          oep(1)=rlread(p(3))
          oep(2)=rlread(p(4))
          oep(3)=rlread(p(5))
          useoep = .true.
          bt_mwl = 0.0d0
          oeptype = 2 !Default to Yukawa potential
          !Read the parameter for the maximum acceptable weight loss
          !for any single eigenstate when truncating the basis set
          if(n.ge.6)bt_mwl = rlread(p(6))
          if(bt_mwl < 0.0d0 .or. bt_mwl > 1.0d0)then
             write(6,*)'Error: BT_MWL invalid in OEP keyword'
             call parallel_end()
          endif
          if(n.ge.7)oeptype = intread(p(7))
          if(n.ge.8)then
             if(intread(p(8))>0)useiterdiag=.true.
          endif
          if(n.ge.9)then
             if(intread(p(9))>0)applyhwithfft=.true.
          endif
!     find the maximum occupied eigenstate for any spin channel
!     Currently not spin polarized
          ie = ifind('e',pname,ntypes)
          do i=1,nspins(ie)
             nocc = max(nocc,nppss(i,ie))
          enddo
! kdelaney: This is like OEP (above) but without the pot options.
!           The potential used for e-p is a bare Coulomb.
       elseif(p(1).eq.'OEPBARE') then ! band potential
          if(ndim.ne.3 .and. root)then
             write(6,*)'ERROR: Findbands currently only for 3D'
             call parallel_end()
          endif
          ecut = rlread(p(2))
          oep(1)=0.0d0
          oep(2)=0.0d0
          oep(3)=0.0d0
          oeptype = 2 !Default to Yukawa potential
          useoep = .true.
          bt_mwl = 0.0d0
          !Read the parameter for the maximum acceptable weight loss
          !for any single eigenstate when truncating the basis set
          if(n.ge.3)bt_mwl = rlread(p(3))
          if(bt_mwl < 0.0d0 .or. bt_mwl > 1.0d0)then
             write(6,*)'Error: BT_MWL invalid in OEP keyword'
             call parallel_end()
          endif
          if(n.ge.4)then
             if(intread(p(4))>0)useiterdiag=.true.
          endif
          if(n.ge.5)then
             if(intread(p(4))>0)applyhwithfft=.true.
          endif
!     find the maximum occupied eigenstate for any spin channel
!     Currently not spin polarized
          ie = ifind('e',pname,ntypes)
          do i=1,nspins(ie)
             nocc = max(nocc,nppss(i,ie))
          enddo
       elseif(p(1).eq.'DFT') then ! Adapt OEP for DFT-LDA bands
          if(.not.useoep)then
             write(6,*)'ERROR: OEP must be used with DFT to gen. 
     .                  trial density. Specify OEP before DFT in
     .                  input file.'
             call parallel_end()
          endif
          !We're using Kerker (see below) for density mixing.
          !A value close to 1 is best for dftmixing in this case.
          dftmixing = rlread(p(2))
          usedft = .true.
          !By default, use Kerker density mixing: best balance of speed
          !and stability.
          scfmethod = 3 !No option to change this for now.
          !Mixmaxg is for Broyden mixing (scfmethod=2). Positive value =>
          !Broyden jacobian is truncated at some G. Unused for all other SCFs
          mixmaxg = -1.0d0
          !If the energy gap (around EF) in the eigenvalue spectrum is small
          !a T=0K Fermi-Dirac occupation will lead to instabilities in the
          !SCF loop. We smear the occupations (which become partial) 
          !with a method that permits subtraction of the entropic contribution
          !to recover T=0 properties. This is only used in finding the 
          !self-consistent DFT density...the Slater determinant is 
          !constructed in the usual way...
          !There are many smearing schemes in findbands, but MP1 is the most
          !reliable.
          smeartype = 2 !Methfessel-Paxton with N=1
          !This is the width of the smearing in Ha. Best values lie between
          !0.005 and 0.04. Too small => Fermi-Dirac at T=0 => oscillations in 
          !SCF again. Too large => many "unoccupied" states required.
          smearwidth = 0.010d0
          if(n>2)then
             smearwidth = rlread(p(3))!Default should be OK in almost all cases
          endif
       elseif(p(1).eq.'JASTROW') then !enable/disable Jastrows
          if(intread(p(2)) .eq. 0) then
             epjastrow = .false.
             cusp_corr = .false.
             if(root)write(6,*)'e-p RPA Jastrow DISABLED'
             if(root)write(6,*)
     .               'OEP cusp correction for finite-ecut: OFF'
          else
             epjastrow = .true.
             cusp_corr = .true.
             if(root)write(6,*)'e-p RPA Jastrow ENABLED'
             if(root)write(6,*)'OEP cusp correction for finite-ecut: ON'
          endif
          if(n.ge.3)then
             if(intread(p(3)) .eq. 0)then
                eejastrow = .false.
                if(root)write(6,*)'e-e RPA Jastrow DISABLED'
             else
                eejastrow = .true.
                if(root)write(6,*)'e-e RPA Jastrow ENABLED'
             endif
          endif
          if(root)then
             write(6,*)
             write(6,*)
          endif
! kdelaney: end changes
       elseif(p(1).eq.'MOLECULES') then
           nmolecules=-1
           molmoves=1
           pc(1)=rlread(p(2))
           pb(1)=1.d0
           if(n.ge.3)bond_length=rlread(p(3))
           if(n.ge.4)pc(2)=rlread(p(4))
           if(n.ge.5)pb(2)=rlread(p(5))
           if(pc(1).gt.0) norbits=1
           if(pc(2).gt.0) norbits=2
           if(norbits.gt.0) then
           if (root) write (6,*)
     &      ' setting up localized functions for electrons ',norbits
           endif
           otype(1)=0 !s state
           otype(2)=1 !pstate
!          open(40,file=qid(1:lpx)//'.bnd'
!    &      ,form='formatted',status='unknown')
       elseif(p(1).eq.'BACKFLOW-ep') then
          ifbackn=intread(p(2))
          if(ifbackn.eq.0) then
             if (root) write(6,*)  ' ep backflow is turned off!'
          else
             if (root) write(6,*)  ' ep backflow is turned on!'
             ifbackflow=1
             plambdab=rlread(p(3))
             psb=rlread(p(4))
             prb=rlread(p(5))
             pwb=rlread(p(6))
             if (root) write(6,*) 'lamb,sb,rb,wb: ',plambdab,psb,prb,pwb
          endif   
!cp
       elseif(p(1).eq.'BACKFLOW-ee') then
          ifbacknee=intread(p(2))
          if(ifbacknee.eq.0) then
             if (root) write(6,*)  ' ee backflow is turned off!'
          else
             if (root) write(6,*)  ' ee backflow is turned on!'
             ifbackflow=1
             plambdabee=rlread(p(3))
             psbee=rlread(p(4))
             prbee=rlread(p(5))
             pwbee=rlread(p(6))
             if (root) 
     &        write(6,*) 'lamb,sb,rb,wb: ',plambdabee,psbee,prbee,pwbee
          endif
       elseif(p(1).eq.'THREEBODY-ep') then
          if3bodyn=intread(p(2))
          if(n.ne.5) then 
             if (root) write(6,*) 'Threebody:# of parameters no match!'
             call parallel_end()
          endif
          if(if3bodyn.eq.0) then
             if (root) write(6,*) ' ep 3-body factor turned off!'
          else
             if (root) write(6,*) ' ep 3-body factor turned on!'
             if3body=1
             plambdat=rlread(p(3))             
             prt=rlread(p(4))
             pwt=rlread(p(5))
             if (root) write(6,*) 'lambT,rT,wT: ',plambdat,prt,pwt
          endif
       elseif(p(1).eq.'THREEBODY-ee') then
          if3bodynee=intread(p(2))
          if(n.ne.5) then
             if (root) write(6,*) 'Threebody:# of parameters no match!'
             call parallel_end()
          endif
          if(if3bodynee.eq.0) then
             if (root) write(6,*) ' ee 3-body factor turned off!'
          else
             if (root) write(6,*) ' ee 3-body factor turned on!'
             if3body=1
             plambdatee=rlread(p(3))
             prtee=rlread(p(4))
             pwtee=rlread(p(5))
             if (root) write(6,*) 'lambT,rT,wT: ',plambdatee,prtee,pwtee
          endif
! mh
       elseif(p(1).eq.'BACKFLOW-A') then
          ifbacka=intread(p(2))
          if(ifbacka.eq.0) then
            if(root) write(6,*)  ' backflow analytical is turned off!'
          else
            if(root) write(6,*)  ' backflow analytical is turned on!'
             ifbackflow=1
             nbact=intread(p(3))
             if(n.ge.4) then
                pla=rlread(p(4))
                if(n.ge.6) then
                  pbee=rlread(p(5))
                  pbep=rlread(p(6))
                end if
             end if
            if(root) write(6,*) 'pla, pbee, pbep ',pla,pbee,pbep
            if(root)write(6,*)'number of k-vec for longrange, nbact='
     .                   ,nbact
          endif
       elseif(p(1).eq.'SELF-A') then
          ifselfa=intread(p(2))
          if(ifselfa.eq.0) then
             if (root) write(6,*)  ' analytical e-e jastrow (self) off'
          else
             if (root) write(6,*)  ' analytical e-e jastrow (self) on '
             if(n.ge.3) then
               pls=rlread(p(3))
             end if
             if (root) write(6,*) 'pls: ',pls
          endif
       elseif(p(1).eq.'EXTRA-A') then
          ifextraa=intread(p(2))
          if(ifextraa.eq.0) then
             if (root) write(6,*)  ' analytical e-p jastrow (extra) off'
          else
             if (root) write(6,*)  ' analytical e-p jastrow (extra) on '
             if(n.ge.3) then
               ple=rlread(p(3))
             end if
             if (root) write(6,*) 'ple: ',ple
          endif
       elseif(p(1).eq.'THREEBODY-A') then
          if3bodya=intread(p(2))
          if(if3bodya.eq.0) then
             if (root) write(6,*)
     .           ' 3-body factor analytical turned off!'
          else
             if (root) write(6,*)
     .          ' 3-body factor analytical turned on!'
             n3act=intread(p(3))
             if(n.ge.4) then
                plambdat=rlread(p(4))
              if(n.ge.6) then
                p3ee=rlread(p(5))
                p3ep=rlread(p(6))
              end if
             end if
             if (root) write(6,*)
     .              'lambT, p3ee, p3ep',plambdat,p3ee,p3ep
             if (root) write(6,*)
     .              'n3act=',n3act
             if3body=1
          endif
! end mh
!      elseif(p(1).eq.'HEDF2') then
!         hedf2=1
!         if(n.ge.2) cutr=rlread(p(2))
       elseif(p(1).eq.'RHO') then
          rho=rlread(p(2))
          if (root) write(6,60) rho
 60       format('Rho  ',f12.6,i4)
       elseif(p(1).eq.'VERBOSE') then
          verbose=intread(p(2))
          if (root) write(6,*) 'verbose=',verbose
       elseif(p(1).eq.'DEBUG') then
          debug=.true.
          idebug=1
          if(root) write(6,*) 'debug turned on!'
       elseif(p(1).eq.'PRESSURE') then
          npt=.true.
          inpt=intread(p(2))
! inpt=1 isotropic change of the box
! inpt=2 rombohedral change (only sides)
! inpt=3 general change (even angles)
          pext=rlread(p(3))
       endif

      go to 1
11    if (root) write (6,*)' finished first pass'
      rewind (1)

! initialize temperature?
!     if(temperature.eq.0.d0) then
!      gState=1
!      beta=1e8
!      if (root) write(6,*) 'Zero temperature calculation, beta set to '
!    +           ,beta
!     else
!      beta=1.d0/temperature
!      gState=0
!      if (root) write(6,*) 'Finite temperature, beta=',beta
!     endif

      if(ntypes.le.0) then
       if (root) write (6,*) 'ntypes.le.0'
       call parallel_end()
      endif
        
!     if(debug) write(70+myrank,'(i5,3g20.10)') 
!    &          (i,(rsc(l,i),l=1,ndim),i=1,nparts)
!     if(debug) write(70+myrank,*)'mctas before setbox: cutr=',cutr(1)
      call setbox(rsc) ! set some box values, k shells
      if(cuta.gt.cutr(1)) then
       if (root) write(6,*) 'cuta.gt.cutr'
       call parallel_end()
      endif
!     if (debug) write(70+myrank,*) 'mctas:after setbox'

      call kset(ndim,kvec,hboxi(1,1,1),volume(1),theta,zerotwist,kmod
     .         ,nkmd,root,debug)
!     if (debug) then
!      write(70+myrank,*) 'mctas:after kset'
!      write(6,*)' mctas: after kset'
!      do is=1,nkact(1)
!       do ik=kmult(is-1,1)+1,kmult(is,1)
!        write(6,'(2i5,3g20.10)') is,ik,(rkcomp(l,ik,1),l=1,ndim)
!       enddo
!      enddo
!      pause
!     endif

      do la=1,ndim
       xold(ltheta+la)=theta(la,1)
      enddo
      xold(ltheta+ndim+1)=phweight(1) 

      if(npslices.gt.1) call addpot(padd,nadd)    !PIMC

      if(rs.gt.0)call coulinit(1) !COUL     
!     if (debug) write(70+myrank,*) 'mctas:after coulinit'

!     if(hedf2.eq.1) call heinit(usr,ulr,lptable,dtable,rho,ndim,cutr
!    .  ,nparts)
!   mh
      iself=2*ifind('e',pname,ntypes)
!   mh end
      if (root) write (6,*) 'mctas:',(p(i),i=1,n)
      call ptrans(0,1,p,1)  ! initialize tables
!     if (debug) write(70+myrank,*) 'mctas:after ptrans'

      lgrad=lx+ndim*nparts
      lcgrad=lgrad+ndim*nparts ! for complex gradient of determinant
      lbuff=lx+4*ndim*nparts  ! size of buffer for configurations
      xold(lage)=0
! put reduced coordinates in stemp for calpps
      ip=ifind('p',pname,ntypes)
      ie=ifind('e',pname,ntypes)
      do i=1,ncomps(ip)
       do l=1,ndim
        stemp(l,i+nfpty(ip)-1)=proton(l,i,1,1)
        stemp(l,i+nfpty(ie)-1)=electron(l,i,1,1)
        do jp=1,npslices
         proton(l,i,jp,2)=proton(l,i,jp,1)
        enddo
       enddo
      enddo
      call calpps(stemp,xold,2,1,1,1) ! initialize configuration properties
      do je=1,neslices
       psit(1,je)=xold(ltf)
       psit(2,je)=xold(ltf)
      enddo
!     if (debug) write(70+myrank,*) 'mctas:after calpps'
      if(root)write (6,*)'initial potential energy ',xold(lv),xold(ltf)
!     if(debug)write(70+myrank,*)
!    &         'initial potential energy ',xold(lv),xold(ltf)
      call initial_av
!      if(root)write(6,*) 'mctas:after initial_av '
      call zeroav(0)
!     if(debug) write(70+myrank,*) 'mctas:after initial_av and zeroav'

! NPT change: 6/6/02
!     if (debug) then
!      write(70+myrank,*) 'box and inverse box matrices'
!      write(70+myrank,'(3f10.4)')((hbox(l,i,1),l=1,ndim),i=1,ndim)
!      write(70+myrank,'(3f10.4)')((hboxi(l,i,1),l=1,ndim),i=1,ndim)
!     endif
      if (jlindem.eq.1) then   ! reduce psites into the unit cube
       it=ifind('p',pname,ntypes)
       do i=1,ncomps(it)
        call reduce(psites(1,i),pbb,dum,hboxi(1,1,1),ndim)
        do l=1,ndim
         psites(l,i)=pbb(l)
        enddo
       enddo
      endif
! initializing array in pt_new address (for constant volume too)
      call setbox_npt(2)
!     if (debug) write(70+myrank,*) 'mctas:after setbox_npt'
      call coulinit(2)
!     if (debug) write(70+myrank,*) 'mctas:after coulinit(2)'

      call run(restart,frestart,qid,timer,time_initial) 

      if (root) write (6,402) 
402   format(' end of simulation')

      call parallel_end()

      end