double precision function chisquare_force(n,x,g,ider)
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
#ifdef MPI
      include 'mpif.h'
      include "COMMON.MPI"
      integer ierror
#endif
      include "COMMON.IOUNITS"
      include "COMMON.FMATCH"
      include "COMMON.CLASSES"
      include "COMMON.WEIGHTS"
      include "COMMON.VMCPAR"
      include "COMMON.PROTFILES"
      include "COMMON.MD"
      include "COMMON.DERIV"
      include "COMMON.INTERACT"
      include "COMMON.ENERGIES"
      include "COMMON.CHAIN"
      include "COMMON.NAMES"
      integer n
      double precision x(n),g(n),gg(max_paropt)
      integer ider
      logical lder
      integer i,ii,iii,ind,jj,j,k,kk,l,ll,iprot,nforce
      double precision aux,chisquare_(maxprot),
     & g_(max_paropt),fmean,fmean_,fsquare,fsquare_,fcov,fcov_
      integer ipass_conf,istart_conf,iend_conf
      double precision tcpu_ini,tcpu
      double precision energia (0:max_ene)
      double precision force_(3,maxres2,maxstr)
#ifdef DEBUG
      write (iout,*) "Enter chisquare_forces, ider",ider
      write (iout,*) "x:",x
      write (iout,*) "ww",ww(:n_ene)
      call flush(iout)
#endif
      calc_grad=.true.
#ifdef MPI
c      write (iout,*) "Before broadcast ider",ider,Master
c      call flush(iout)
      call MPI_Bcast(ider,1,MPI_INTEGER,Master,Comm1,ierror)
c      write (iout,*) "After broadcast ider ",ider
c      call flush(iout)
#endif
      if (ider.eq.0) return
      lder=ider.eq.2

      if (lder) g=0.0d0

      do iprot=1,nprot

      chisquare(iprot)=0.0d0
      
      gg = 0.0d0

      fmean=0.0d0
      fsquare=0.0d0
      fcov=0.0d0

      if (.not.fmatch(iprot)) cycle

      call restore_molinfo(iprot)

      if (.not. mod_other_params) then
        open (ientin,file=bforcefiles(iprot),status="old",
     &    form="unformatted",access="direct",recl=lenrec_forces(iprot))
#ifdef DEBUG
        write (iout,*) "bforcefiles",bforcefiles(iprot)
#endif
      else
        open (icbase,file=bprotfiles_MD(iprot),status="old",
     &    form="unformatted",access="direct",recl=lenrec_MD(iprot))
      endif
c      write (iout,*) "chisquare: files opened"
c      call flush(iout)
      ipass_conf=0
#ifdef MPI
      do i=indstart_MD(me1,iprot),indend_MD(me1,iprot),maxstr_proc
        istart_conf=i
        iend_conf=min0(i+maxstr_proc-1,indend_MD(me1,iprot))
#else
      do i=1,ntot_work_MD(iprot),maxstr_proc
        istart_conf=i
        iend_conf=min0(i+maxstr_proc-1,ntot_MD(iprot))
#endif
c        write (iout,*) "i",i," istart_conf",istart_conf," iend_conf",
c     &    iend_conf
c        call flush(iout)
        ii=0
c
c Read the chunk of conformations off a DA scratchfile.
c
        if (.not. mod_other_params) then
#ifdef DEBUG
          write (iout,*) "chisquare: calling daread_forces",istart_conf,
     &      iend_conf
#endif
          call daread_forces(iprot,istart_conf,iend_conf)
          do iii=istart_conf,iend_conf
            ii=ii+1
            do k=1,nsite(iprot)
              do l=1,3
                aux=0.0d0
                do j=1,n_ene
                  aux=aux+forctb(j,l,k,ii,iprot)*ww(j)
                enddo
                force(l,k,iii,iprot)=aux
              enddo
            enddo
          enddo
        else
          call daread_MDcoords(iprot,istart_conf,iend_conf)
          do iii=istart_conf,iend_conf
            ii=ii+1
            call restore_MDcoords(iprot,ii)
            call int_from_cart1(.false.)
#ifdef DEBUG
            write (iout,*) "Force: CG Coordinates conf",iii,ii
            do k=1,nres
              if (itype(k).eq.10 .or. itype(k).eq.ntyp1) then
                write (iout,'(a4,i5,3f10.3)')
     &            restyp(itype(k)),k,(c(j,k),j=1,3)
                else
                  write (iout,'(a4,i5,3f10.3,5x,3f10.3)')
     &            restyp(itype(k)),k,(c(j,k),j=1,3),(c(j,k+nres),j=1,3)
               endif
            enddo
            call int_from_cart1(.true.)
#endif
            call zerograd
            call etotal(energia(0))
            call cartgrad
            do k=1,nres
              do l=1,3
                force(l,k,iii,iprot)=-gcart(l,k)
                do j=1,n_ene
                  forctb(j,l,k,ii,iprot)=-gcompon(j,l,k)
                enddo
              enddo
            enddo
            kk=nres
            do k=nnt,nct
              if (itype(k).eq.10 .or. itype(k).eq.ntyp1) cycle
              kk=kk+1
              do l=1,3
                force(l,kk,iii,iprot)=-gxcart(l,k)
                do j=1,n_ene
                  forctb(j,l,kk,ii,iprot)=-gcomponx(j,l,k)
                enddo
              enddo
            enddo 
#ifdef CHECKFORCE
            write (iout,*) "Checking forces from components"
            write (iout,'(2a5,3a10)')"coor","site","full",
     &          "from compon","diff"
            do k=1,nsite(iprot)
              do l=1,3
                aux=0.0d0
                do j=1,n_ene
                  aux=aux+forctb(j,l,k,ii,iprot)*ww(j)
                enddo
                write(iout,'(2i5,3f10.5)') l,k,force(l,k,iii,iprot),aux,
     &             (force(l,k,iii,iprot)-aux)/force(l,k,iii,iprot)
              enddo
            enddo
#endif
          enddo
        endif
c Compute chisquae components
        ii=0
        do iii=istart_conf,iend_conf
#ifdef DEBUG
          write (iout,*) "Forces: conformation",iii,ii
          jj=nres
          do j=1,nres
            if (itype(j).eq.10 .or. itype(j).eq.ntyp1) then
              write (iout,'(a4,i5,3f10.5)') 
     &          restyp(itype(j)),j,(force(k,j,iii,iprot),k=1,3)
            else
              jj=jj+1
              write (iout,'(a4,i5,3f10.5,5x,3f10.5)')
     &          restyp(itype(j)),j,(force(k,j,iii,iprot),k=1,3),
     &           (force(k,jj,iii,iprot),k=1,3)
            endif
          enddo
#endif
          ii=ii+1
          do k=1,nsite(iprot)
            do l=1,3
              chisquare(iprot)=chisquare(iprot)+
     &          (CGforce_MD(l,k,iii,iprot)-force(l,k,iii,iprot))**2
              fmean=fmean+force(l,k,iii,iprot)
              fsquare=fsquare+force(l,k,iii,iprot)**2
              fcov=fcov+force(l,k,iii,iprot)*CGforce_MD(l,k,iii,iprot)
            enddo
          enddo
          if (lder) then
          ll=0
          do l=1,n_ene
            if (imask(l).gt.0) then
            ll=ll+1
            do j=1,nsite(iprot)
              do k=1,3
                gg(ll)=gg(ll)+(CGforce_MD(k,j,iii,iprot)
     &            -force(k,j,iii,iprot))*forctb(l,k,j,ii,iprot)
c                write (iout,*) "iii",iii," ii",ii,"l",l," ll",ll," j",j
c                write (iout,*) " residual",
c     &           CGforce_MD(k,j,iii,iprot)-force(k,j,iii,iprot),
c     &           " forctab",forctb(l,k,j,ii,iprot)
c                write (iout,*) "gg",gg(ll)
              enddo
            enddo  
            endif
          enddo
          endif
        enddo
      enddo

      if (lder) then
        ll=0
        do l=1,n_ene
         if (imask(l).gt.0) then
           ll=ll+1
           g(ll)=g(ll)+wforce(iprot)*gg(ll)/
     &       (3*ntot_work_MD(iprot)*nsite(iprot))
         endif
        enddo
c        write (iout,*) "chisquare ll:",ll
      endif
c      write (iout,*) "fmean",fmean," fsquare",fsquare
#ifdef MPI
      fmean_=fmean
      call MPI_Reduce(fmean_,fmean,1,MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,Comm1,ierror)
      fsquare_=fsquare
      call MPI_Reduce(fsquare_,fsquare,1,MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,Comm1,ierror)
      fcov_=fcov
      call MPI_Reduce(fcov_,fcov,1,MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,Comm1,ierror)
#endif
      nforce = 3*ntot_work_MD(iprot)*nsite(iprot)
      fmean=fmean/nforce
c      write (iout,*) "ntot_work",ntot_work_MD(iprot)," nforce",nforce
c      write (iout,*) "fmean",fmean/nforce," fsquare",fsquare/nforce,
c     & " variance",(fsquare-nforce*fmean)/nforce
      rcorr(iprot)=(fcov-nforce*meanforce(iprot)*fmean)
     &  /dsqrt(nforce*(fsquare-nforce*fmean*fmean)*varforce(iprot))
      enddo !v iprot
#ifdef DEBUG
      do iprot=1,nprot
        write (iout,*) "Protein",iprot," forces"
        write (iout,"(4x,a,t55,a)") "Backbone","Sidechain"
#ifdef MPI
        do iii=indstart_MD(me1,iprot),indend_MD(me1,iprot)
#else
        do iii=1,ntot_work_MD(iprot)
#endif
          write (iout,*) "Conformation",iii
          ii=nres
          do i=1,nres
            if (itype(i).eq.10 .or.itype(i).eq.ntyp1) then
              write (iout,'(a4,i5,3e15.5)') restyp(itype(i)),i,
     &          (CGforce_MD(j,i,iii,iprot),j=1,3)
              write (iout,'(9x,3e15.5)') (force(j,i,iii,iprot),j=1,3)
            else
              ii=ii+1
              write (iout,'(a4,i5,3e15.5,5x,3e15.5)')restyp(itype(i)),i,
     &          (CGforce_MD(j,i,iii,iprot),j=1,3),
     &          (CGforce_MD(j,ii,iii,iprot),j=1,3)
              write (iout,'(9x,3e15.5,5x,3e15.5)')
     &          (force(j,i,iii,iprot),j=1,3),
     &          (force(j,ii,iii,iprot),j=1,3)
            endif
          enddo
        enddo
      enddo

      write(iout,*)"chisquare",chisquare(:nprot)
#endif
c Compute total chisquare and gradient
#ifdef MPI
      chisquare_=chisquare
      call MPI_Reduce(chisquare_,chisquare,nprot,MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,Comm1,ierror)
      if (lder) then
        gg=g
c        write (iout,*) "g before reduce",g(:n)
        call MPI_Reduce(gg,g,n,MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,Comm1,ierror)
c        write (iout,*) "g after reduce",g(:n)
      endif
#endif
c      write (iout,*) "chisquare",chisquare(:nprot)
      aux=0.0d0
      do iprot=1,nprot
        if (.not.fmatch(iprot)) cycle
        chisquare(iprot)=chisquare(iprot)/
     &    (3*ntot_work_MD(iprot)*nsite(iprot))
        aux=aux+0.5d0*wforce(iprot)*chisquare(iprot)
      enddo
#ifdef DEBUG
      write (iout,*) "Leave chisquare_forces"    
      call flush(iout)
#endif
#ifdef MPI
      if (ider.eq.3) then
      do iprot=1,nprot
#ifdef DEBUG
      write (iout,*) "fforce gather me1",me1," indstart",
     & indstart_MD(me1,iprot)," indend",indend_MD(me1,iprot)
#endif
      do i=indstart_MD(me1,iprot),indend_MD(me1,iprot)
        do j=1,nsite(iprot) 
          do k=1,3
            force_(k,j,i-indstart_MD(me1,iprot)+1)=force(k,j,i,iprot)
            force(k,j,i,iprot)=0.0d0
          enddo
        enddo
      enddo
#ifdef DEBUG
      write (iout,*) "force_ scount_MD",scount_MD(me1,iprot)
      do i=1,scount_MD(me1,iprot)
        do j=1,nsite(iprot)
          write (iout,'(2i5,3f10.5)') i,j,(force_(k,j,i),k=1,3)
        enddo
      enddo
      write (iout,*) "scount_MD",scount_MD(:nprocs-1,iprot)
      write (iout,*) "idispl_MD",idispl_MD(:nprocs-1,iprot)
#endif
      call MPI_Gatherv(force_(1,1,1),scount_MD(me1,iprot),MPI_FORCE,
     &   force(1,1,1,iprot),scount_MD(0,iprot),idispl_MD(0,iprot),
     &   MPI_FORCE,Master,Comm1, IERROR)
#ifdef DEBUG
      write (iout,*) "Gather ierror",ierror
      write (iout,*) "force after gather"
      do i=1,ntot_MD(iprot)
        do j=1,nsite(iprot)
          write (iout,'(2i5,3f10.5)') i,j,(force(k,j,i,iprot),k=1,3)
        enddo
      enddo
#endif
      enddo
      endif ! iderp.eq.3
#endif
      chisquare_force=aux
      return
      end