double precision function maxlik_pdb(n,x,g,ider)
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      integer n,nmax
      parameter (nmax=max_paropt)
#ifdef MPI
      include 'mpif.h'
      include 'COMMON.MPI'
      integer ierror
      double precision g_(max_paropt),sumlik_(0:ntyp,0:ntyp),
     &  partf_,emin_,der_partf_(maxvar)
#endif
      include "COMMON.IOUNITS"
      include "COMMON.PMFDERIV"
      include "COMMON.WEIGHTS"
      include "COMMON.TORSION"
      include "COMMON.FFIELD"
      include "COMMON.PDBSTAT"
      include "COMMON.GEO"
      double precision b1m(3,2),b2m(3,2),dm(2,2),cm(2,2),
     &  b1m2(3,2),b2m1(3,2),b11,b12,b21,b22,c11,c12,d11,d12,c1(0:3),
     &  c2(0:3),em(2,2,2),em2(2,2,2),e0m(3),e0m2(3),jac(nmax)
      double precision v1_kcc_loc(3,3,2),v2_kcc_loc(3,3,2),
     & etoriv1(3,3,2),etoriv2(3,3,2),enei,boltz,emin,
     & energia(36*36*72)
      double precision theta1,theta2,gam,thetap1,thetap2,gam1,gam2,
     & phi,cost1,cost2,sint1,sint2,
     & cosg,sing,cos2g,sin2g,sum,d3,delta,f,fb11,fb21,fb12,fb22,fc11,
     & fc12,fd11,fd12,etori,sint1sint2,cosphi,sinphi,sumvalc,sumvals,
     & vol
      double precision conv /0.01745329251994329576d0/
      double precision r0pp(3)  /5.2739d0,5.4561d0,5.2261d0/
      double precision epspp(3) /1.6027d0,1.4879d0,0.0779d0/
      double precision dCACA(2) /3.784904d0,3.810180d0/
      double precision facp1(4),facpp(4),facturn3,ddist
      logical lder
      double precision x(n),g(n)
      integer ipoint,iipoint,ip,it1,it2,itp1,itp2,i,j,k,l,ii,jj!,irt1,irt2
      integer ider
      double precision sumlik(0:ntyp,0:ntyp),partf,der_partf(maxvar)
      integer itheta1,itheta2,igam
      vol = dlog((incr_theta*deg2rad)**2*incr_gam*deg2rad)
      delta = pi
#ifdef MPI
c nfun.eq.-1 : no broadcast
c n.eq.-1 : stop function evaluation
c nfun.ne.-1 .and. n.ne.-1 : master distribute variables to slaves, all do their share
c of computations
      call MPI_Bcast(ider,1,MPI_LOGICAL,Master,MPI_COMM_WORLD,ierror)
      if (ider.le.0) return
#endif
      lder=ider.gt.1
c      write (iout,*) "MAXLIKPDB: ider",ider," lder",lder
c      call flush(iout)
      IF (lder) THEN
c Zero out gradient and hessian
      do i=1,n
        g(i)=0.0d0
      enddo
      ENDIF

      sum=0.0d0

      sumlik=0.0d0
      do i=1,npoint_pdb
#ifdef DEBUG
      write (iout,*) "point",i," it1",it1," it2",it2
#endif
      it1 = ipdbtyp(1,i)
      it2 = ipdbtyp(2,i)
      theta1=zz_pdb(1,i)
      theta2=zz_pdb(2,i)
      gam=(zz_pdb(3,i)+delta-pi)
      call torval_e_der(n,it1,it2,theta1,theta2,gam,enei,jac,lder)
#ifdef DEBUG
      write (2,*) "i",i," it1",it1," it2",it2," theta1",theta1,
     &  " theta2",theta2," gamma",gam," enei",enei
#endif
      sumlik(it1,it2)=sumlik(it1,it2)+beta_pdb*enei*resol(i)
      IF (lder) THEN
c Contributions to the gradient and the hessian
      do k=1,n
        g(k)=g(k)+beta_pdb*jac(k)*resol(i)
      enddo
      if (mask_beta_PDB.gt.0) 
     &  g(mask_beta_PDB)=g(mask_beta_PDB)+enei*resol(i)
#ifdef DEBUG
      print *,"i",i," jac",jac(:n)," resol",resol(i)," g",g(:n)
#endif
      ENDIF ! lder
      enddo  ! i
#ifdef MPI
      sumlik_=sumlik
      call MPI_Reduce(sumlik_,sumlik,(ntyp+1)*(ntortyp+1),
     &      MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,ierror)
      if (lder) then
        g_=g
        call MPI_Reduce(g_(1),g(1),n,MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,Comm1,ierror)
c        write (iout,*) "MAXLIKPDB: after reducing gradient"
c        call flush(iout)
      endif
#endif
#ifdef DEBUG
      do it1=0,ntortyp-1
        do it2=0,ntortyp-1
          write (iout,*) "sumlik",it1,it2,sumlik(it1,it2)
        enddo
      enddo
#endif
c      sumlik=0
c      if (lder) g=0
#define CIPISZCZE
#ifdef CIPISZCZE
c Compute the partition functions for consecutive pairs of residue types
      do it1=0,ntortyp-1
        do it2=0,ntortyp-1
#ifdef ENEOUT
          write (iout,'(a,2i2,2i10)') 
     &      "BEGIN",it1,it2,istart_calka,iend_calka
#endif
          ii=0
          do itheta1=70,160-incr_theta,incr_theta
            do itheta2=70,160-incr_theta,incr_theta
              do igam=-180,180-incr_gam,incr_gam
                ii=ii+1
                if (ii.ge.istart_calka .and. ii.le.iend_calka) then
                  theta1=(itheta1+0.5d0)*deg2rad
                  theta2=(itheta2+0.5d0)*deg2rad
                  gam=(igam+0.5d0)*deg2rad
                  call torval_e_der(n,it1,it2,theta1,theta2,gam,enei,
     &              jac,.false.)
                  energia(ii)=enei
#ifdef ENEOUT
c                  write (iout,'(3i5,1pe15.5)') itheta1,itheta2,igam,enei
#endif
                endif
              enddo
            enddo
          enddo
          emin=energia(istart_calka)
          do ii=istart_calka+1,iend_calka
            if (energia(ii).lt.emin) emin=energia(ii)
          enddo
#ifdef MPI
          emin_=emin
          call MPI_AllReduce(emin_,emin,1,
     &         MPI_DOUBLE_PRECISION,MPI_MIN,Comm1,ierror)
#endif
c          emin=0.0d0
c          write (2,*) it1,it2," emin",emin
          partf=0.0d0
          der_partf=0.0d0
          ii=0
          do itheta1=70,160-incr_theta,incr_theta
            do itheta2=70,160-incr_theta,incr_theta
              do igam=-180,180-incr_gam,incr_gam
                ii=ii+1
                if (ii.ge.istart_calka .and. ii.le.iend_calka) then
                  theta1=(itheta1+0.5d0)*deg2rad
                  theta2=(itheta2+0.5d0)*deg2rad
                  gam=(igam+0.5d0)*deg2rad
c                    call torval_e_der(n,it1,it2,theta1,theta2,gam,enei,
c     &              jac,lder)
                  enei=energia(ii)-emin
                  boltz = dsin(theta1)*dsin(theta2)*dexp(-beta_pdb*enei)
                  partf=partf+boltz
c                  write (iout,*) "it1",it1," it2",it2,
c     &              " theta1",theta1," theta2",theta2,
c     &              " gamma",gam," boltz",boltz," enei",enei,
c     &              " partf",partf
                  if (lder) then
                    call torval_e_der(n,it1,it2,theta1,theta2,gam,enei,
     &              jac,lder)
                    do i=1,n
                      der_partf(i)=der_partf(i)-boltz*jac(i)
                    enddo
                    if (mask_beta_PDB.gt.0)
     &              der_partf(mask_beta_PDB)=der_partf(mask_beta_PDB)-
     &               boltz*enei/beta_pdb
                  endif
                endif
              enddo
            enddo
          enddo
#ifdef ENEOUT
          write (iout,'(a,2i2)') "END",it1,it2
#endif
c          write (iout,*) "partf before REDUCE",it1,it2,partf
#ifdef MPI
          partf_=partf
          call MPI_Reduce(partf_,partf,1,MPI_DOUBLE_PRECISION,
     &           MPI_SUM,Master,Comm1,ierror)
c          write (iout,*) "Afte REDUCE: partf_",partf_," partf",partf
          if (lder) then
            der_partf_=der_partf
            call MPI_Reduce(der_partf_(1),der_partf(1),n,
     &         MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,ierror)
c            write (iout,*) "MAXLIKPDB: after reducing gradient"
c            call flush(iout)
          endif
          if (me.eq.Master) then
#endif
          sumlik(it1,it2)=sumlik(it1,it2)+dlog(partf)-beta_pdb*emin+vol
#ifdef DEBUG
          write(iout,*) "it1",it1," it2",it2," partf",partf," sumlik",
     &      sumlik(it1,it2)," emin",emin
#endif
          if (lder) then
            do i=1,n
              g(i)=g(i)+beta_pdb*der_partf(i)/partf
            enddo
          endif
#ifdef MPI
        endif
#endif
        enddo
      enddo
#ifdef ENEOUT
      stop
#endif
c      write (iout,*) "MAXLIKPDB: end of loop"
c      call flush(iout)
#endif
c Compute the total contribution to the target function from a given core
      sum=0.0d0
      do it1=0,ntortyp-1
        do it2=0,ntortyp-1
          sum=sum+sumlik(it1,it2)
        enddo
      enddo
      maxlik_pdb = sum
      return
      end
c----------------------------------------------------------------------
      subroutine torval_e_der(nvar,it1,it2,theta1,theta2,gam,enei,jac,
     &  lder)
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      integer n,nmax
      parameter (nmax=max_paropt)
#ifdef MPI
      include 'mpif.h'
      include 'COMMON.MPI'
      integer ierror
      double precision g_(max_paropt),sum_
#endif
      include "COMMON.IOUNITS"
      include "COMMON.PMFDERIV"
      include "COMMON.WEIGHTS"
      include "COMMON.TORSION"
      include "COMMON.FFIELD"
      include "COMMON.GEO"
      integer nvar
      integer it1,it2,itp1,itp2,ip,i,j,k,l,ii,jj
      double precision b1m(3,2),b2m(3,2),dm(2,2),cm(2,2),
     &  b1m2(3,2),b2m1(3,2),b11,b12,b21,b22,c11,c12,d11,d12,c1(0:3),
     &  c2(0:3),em(2,2,2),em2(2,2,2),e0m(3),e0m2(3),jac(nmax)
      double precision v1_kcc_loc(3,3,2),v2_kcc_loc(3,3,2),
     & etoriv1(3,3,2),etoriv2(3,3,2),enei,etori,ebendi1,ebendi2
      double precision theta1,theta2,gam,thetap1,thetap2,gam1,gam2,
     & phi,cost1,cost2,sint1,sint2,
     & cosg,sing,cos2g,sin2g,sum,d3,delta,f,fb11,fb21,fb12,fb22,fc11,
     & fc12,fd11,fd12,sint1sint2,cosphi,sinphi,sumvalc,sumvals
      integer ind_ang /11/, ind_tor /13/
      logical lder
      double precision tschebyshev
c following to test with selected types
      if (it1.le.1) then
        itp1=1
      else
        itp1=2
      endif
      if (iabs(it2).le.1) then
        itp2=1
      else
        itp2=2
      endif

      ip=ip+1

      delta=pi

      do i=1,2
        do k=1,3
          b1m(k,i)=bnew1tor(k,i,it2)
          b1m2(k,i)=bnew2tor(k,i,it2)
          b2m1(k,i)=bnew1tor(k,i,it1)
          b2m(k,i)=bnew2tor(k,i,it1)
        enddo
      enddo
      do i=1,2
        do k=1,2
          cm(k,i)=ccnewtor(k,i,it2)
          dm(k,i)=ddnewtor(k,i,it1)
        enddo
      enddo
      do i=1,2
        do j=1,2
          do k=1,2
            em(k,j,i)=eenewtor(k,j,i,it1)
            em2(k,j,i)=eenewtor(k,j,i,it2)
          enddo
        enddo
      enddo
c Invert coefficients for L-chirality
#ifdef DEBUG
      write(iout,*) "it1",it1," it2",it2
      write(iout,*) "b1i (b1m)"
      do i=1,3
        write(iout,'(2f10.5,5x,2f10.5)') b1m(i,1),b1m(i,2)
      enddo
      write(iout,*) "b2j (b1m2)"
      do i=1,3
        write(iout,'(2f10.5,5x,2f10.5)') b1m2(i,1),b1m2(i,2)
      enddo
      write(iout,*) "b1j (b2m1)"
      do i=1,3
        write(iout,'(2f10.5,5x,2f10.5)') b2m1(i,1),b2m1(i,2)
      enddo
      write(iout,*) "b2j (b2m)"
      do i=1,3
        write(iout,'(2f10.5,5x,2f10.5)') b2m(i,1),b2m(i,2)
      enddo
      write (iout,*) "c"
      do i=1,2
        write (iout,'(2f10.5)') (cm(i,j),j=1,2)
      enddo
      write (iout,*) "d"
      do i=1,2
        write (iout,'(2f10.5)') (dm(i,j),j=1,2)
      enddo
c      print *,"em1",em," em2",em2," e0m",e0m," e02m",e0m2
#endif
      v1_kcc_loc=0.0d0
      v2_kcc_loc=0.0d0
      etoriv1=0.0d0
      etoriv2=0.0d0
      do k=1,3
        do l=1,3
          v1_kcc_loc(k,l,1)=b2m(l,1)*b1m(k,1)+b2m(l,2)*b1m(k,2)
          v2_kcc_loc(k,l,1)=b2m(l,2)*b1m(k,1)+b2m(l,1)*b1m(k,2)
        enddo
      enddo
      do k=1,2
        do l=1,2
          v1_kcc_loc(k,l,2)=-(cm(k,1)*dm(l,1)-cm(k,2)*dm(l,2))
          v2_kcc_loc(k,l,2)=-(cm(k,2)*dm(l,1)+cm(k,1)*dm(l,2))
        enddo
      enddo
      cost1=dcos(theta1)
      cost2=dcos(theta2)
      sint1=dsin(theta1)
      sint2=dsin(theta2)
      cosg=dcos(gam)
      sing=dsin(gam)
      cos2g=dcos(2*gam)
      sin2g=dsin(2*gam)
c check with the formula from unres
      c1(0)=0.0d0
      c2(0)=0.0d0
      c1(1)=1.0d0
      c2(1)=1.0d0
      do j=2,3
        c1(j)=c1(j-1)*cost1
        c2(j)=c2(j-1)*cost2
      enddo
      etori=0.0d0
      sint1sint2=1.0d0
      do j=1,2
        cosphi=dcos(j*gam)
        sinphi=dsin(j*gam)
        sint1sint2=sint1sint2*sint1*sint2
        sumvalc=0.0d0
        sumvals=0.0d0
        do k=1,3
          do l=1,3
            sumvalc=sumvalc+v1_kcc_loc(l,k,j)*c1(k)*c2(l)
            sumvals=sumvals+v2_kcc_loc(l,k,j)*c1(k)*c2(l)
            etoriv1(l,k,j)=c1(k)*c2(l)*sint1sint2*cosphi
            etoriv2(l,k,j)=c1(k)*c2(l)*sint1sint2*sinphi
          enddo
        enddo
        etori=etori+(sumvalc*cosphi+sumvals*sinphi)*sint1sint2
      enddo
c Bending contributions
      ebendi1=v1bend_chyb(0,it1)+
     &    tschebyshev(1,nbend_kcc_Tb(it1),v1bend_chyb(1,it1),cost1)
      ebendi2=v1bend_chyb(0,it2)+
     &    tschebyshev(1,nbend_kcc_Tb(it2),v1bend_chyb(1,it2),cost2)
c      write (iout,*) "etori",etori," ebendi1",ebendi1," ebendi2",ebendi2
c Compute derivatives of the energy in torsional parameters
      enei = wtor*etori+wang*(ebendi1+ebendi2)
#ifdef ENEOUT
      write (iout,'(3f8.1,4(1pe15.5))') theta1*rad2deg,theta2*rad2deg,
     &    gam*rad2deg,etori,ebendi1,ebendi2,enei
#endif
c
      IF (lder) THEN
c
      do k=1,nvar
        jac(k)=0.0d0
      enddo
c Free term
c Derivatives in b1 of the second residue
c        jj=-1
      do j=1,3
c          jj=jj+2
        jj=iabs(mask_bnew1tor(j,1,iabs(it2)))
        if (jj.gt.0) then
          do k=1,3
            jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,1)
     &        +etoriv2(j,k,1)*b2m(k,2)
          enddo
        endif
        jj=iabs(mask_bnew1tor(j,2,iabs(it2)))
        if (jj.gt.0) then
          do k=1,3
            if (it2.gt.0) then
              jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,2)
     &        +etoriv2(j,k,1)*b2m(k,1)
            else if (it2.lt.0) then
              jac(jj)=jac(jj)-etoriv1(j,k,1)*b2m(k,2)
     &        -etoriv2(j,k,1)*b2m(k,1)
            endif
          enddo
        endif
      enddo
c Derivatives in b2 of the first residue
      do j=1,3
c          jj=jj+2
        jj=iabs(mask_bnew2tor(j,1,it1))
        if (jj.gt.0) then
          do k=1,3
            jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,1)
     &      +etoriv2(k,j,1)*b1m(k,2)
          enddo
        endif
        jj=iabs(mask_bnew2tor(j,2,it1))
        if (jj.gt.0) then
          do k=1,3
            if (it1.gt.0) then
              jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,2)
     &        +etoriv2(k,j,1)*b1m(k,1)
            endif
          enddo
        endif
      enddo
c Derivatives in c
      do j=1,2
c        jj=jj+2
        jj=iabs(mask_ccnewtor(j,1,iabs(it2)))
        if (jj.gt.0) then
          do k=1,2
            jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,1)
     &       -etoriv2(j,k,2)*dm(k,2)
          enddo
        endif
        jj=iabs(mask_ccnewtor(j,2,iabs(it2)))
        if (jj.gt.0) then
          do k=1,2
            if (it2.gt.0) then
              jac(jj)=jac(jj)+etoriv1(j,k,2)*dm(k,2)
     &        -etoriv2(j,k,2)*dm(k,1)
            else if (it2.lt.0) then
              jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,2)
     &        +etoriv2(j,k,2)*dm(k,1)
            endif
          enddo
        endif
      enddo
c Derivatives in d
      do j=1,2
        jj=iabs(mask_ddnewtor(j,1,it1))
        if (jj.gt.0) then
          do k=1,2
            jac(jj)=jac(jj)-etoriv1(k,j,2)*cm(k,1)
     &       -etoriv2(k,j,2)*cm(k,2) 
          enddo
        endif
        jj=iabs(mask_ddnewtor(j,2,it1))
        if (jj.gt.0) then
          do k=1,2
            if (it1.gt.0) then
              jac(jj)=jac(jj)+etoriv1(k,j,2)*cm(k,2)
     &        -etoriv2(k,j,2)*cm(k,1)
            endif
          enddo
        endif
      enddo
      do i=1,nvar
        jac(i)=wtor*jac(i)
      enddo
c      jac=0.0d0
      do j=1,nbend_kcc_TB(it1)
        ii = mask_v1bend_chyb(j,it1)
        if (ii.gt.0) jac(ii) = wang*dcos(j*theta1)
      enddo
      do j=1,nbend_kcc_TB(it2)
        ii = mask_v1bend_chyb(j,it2)
        if (ii.gt.0) jac(ii) = jac(ii) + wang*dcos(j*theta2)
      enddo
      if (imask(ind_ang).gt.0) jac(imask(ind_ang))=ebendi1+ebendi2
      if (imask(ind_tor).gt.0) jac(imask(ind_tor))=etori

      ENDIF

      return
      end