subroutine initialize
C 
C Define constants and zero out tables.
C
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'sizesclu.dat'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.TORSION'
      include 'COMMON.FFIELD'
      include 'COMMON.SBRIDGE'
      include 'COMMON.MINIM' 
      include 'COMMON.DERIV'
      include "COMMON.NAMES"
      include "COMMON.TIME1"
C
C The following is just to define auxiliary variables used in angle conversion
C
      pi=4.0D0*datan(1.0D0)
      dwapi=2.0D0*pi
      dwapi3=dwapi/3.0D0
      pipol=0.5D0*pi
      deg2rad=pi/180.0D0
      rad2deg=1.0D0/deg2rad
      angmin=10.0D0*deg2rad
      Rgas = 1.987D-3
C
C Define I/O units.
C
      inp=    1
      iout=   2
      ipdbin= 3
      ipdb=   7
      imol2= 18
      jplot= 19
      jstatin=10
      imol2=  4
      igeom=  8
      intin=  9
      ithep= 11
      irotam=12
      itorp= 13
      itordp= 23
      ielep= 14
      isidep=15 
      isidep1=22
      iscpp=25
      icbase=16
      ifourier=20
      istat= 17
      ibond=28
      isccor=29
      jrms=30
C
C Set default weights of the energy terms.
C
      wlong=1.0D0
      welec=1.0D0
      wtor =1.0D0
      wang =1.0D0
      wscloc=1.0D0
      wstrain=1.0D0
C
C Zero out tables.
C
      ndih_constr=0
      do i=1,maxres2
	do j=1,3
	  c(j,i)=0.0D0
	  dc(j,i)=0.0D0
        enddo
      enddo
      do i=1,maxres
	do j=1,3
	  xloc(j,i)=0.0D0
        enddo
      enddo
      do i=1,ntyp
	do j=1,ntyp
	  aa(i,j)=0.0D0
	  bb(i,j)=0.0D0
	  augm(i,j)=0.0D0
	  sigma(i,j)=0.0D0
	  r0(i,j)=0.0D0
	  chi(i,j)=0.0D0
        enddo
	do j=1,2
	  bad(i,j)=0.0D0
        enddo
	chip(i)=0.0D0
	alp(i)=0.0D0
	sigma0(i)=0.0D0
	sigii(i)=0.0D0
	rr0(i)=0.0D0
	a0thet(i)=0.0D0
	do j=1,2
         do ichir1=-1,1
          do ichir2=-1,1
          athet(j,i,ichir1,ichir2)=0.0D0
          bthet(j,i,ichir1,ichir2)=0.0D0
          enddo
         enddo
        enddo
	do j=0,3
	  polthet(j,i)=0.0D0
        enddo
	do j=1,3
	  gthet(j,i)=0.0D0
        enddo
	theta0(i)=0.0D0
	sig0(i)=0.0D0
	sigc0(i)=0.0D0
	do j=1,maxlob
	  bsc(j,i)=0.0D0
	  do k=1,3
	    censc(k,j,i)=0.0D0
          enddo
          do k=1,3
	    do l=1,3
	      gaussc(l,k,j,i)=0.0D0
            enddo
          enddo
	  nlob(i)=0
        enddo
      enddo
      nlob(ntyp1)=0
      dsc(ntyp1)=0.0D0
      do i=-maxtor,maxtor
	itortyp(i)=0
       do iblock=1,2
        do j=-maxtor,maxtor
          do k=1,maxterm
            v1(k,j,i,iblock)=0.0D0
            v2(k,j,i,iblock)=0.0D0
           enddo
         enddo      
       enddo
      enddo
      do iblock=1,2
       do i=-maxtor,maxtor
        do j=-maxtor,maxtor
         do k=-maxtor,maxtor
          do l=1,maxtermd_1
            v1c(1,l,i,j,k,iblock)=0.0D0
            v1s(1,l,i,j,k,iblock)=0.0D0
            v1c(2,l,i,j,k,iblock)=0.0D0
            v1s(2,l,i,j,k,iblock)=0.0D0
          enddo !l
          do l=1,maxtermd_2
           do m=1,maxtermd_2
            v2c(m,l,i,j,k,iblock)=0.0D0
            v2s(m,l,i,j,k,iblock)=0.0D0
           enddo !m
          enddo !l
        enddo !k
       enddo !j
      enddo !i
      enddo !iblock
      do i=1,maxres
	itype(i)=0
	itel(i)=0
      enddo
C Initialize the bridge arrays
      ns=0
      nss=0 
      nhpb=0
      do i=1,maxss
	iss(i)=0
      enddo
      do i=1,maxss
	dhpb(i)=0.0D0
      enddo
      do i=1,maxss
	ihpb(i)=0
	jhpb(i)=0
      enddo
C
C Initialize timing.
C
      call set_timers
C
C Initialize variables used in minimization.
C   
c     maxfun=5000
c     maxit=2000
      maxfun=500
      maxit=200
      tolf=1.0D-2
      rtolf=5.0D-4
C 
C Initialize the variables responsible for the mode of gradient storage.
C
      nfl=0
      icg=1
      do i=1,14
        do j=1,14
          if (print_order(i).eq.j) then
            iw(print_order(i))=j
            goto 1121
          endif
        enddo
1121    continue
      enddo
      calc_grad=.false.
C Set timers and counters for the respective routines
      t_func = 0.0d0
      t_grad = 0.0d0
      t_fhel = 0.0d0
      t_fbet = 0.0d0
      t_ghel = 0.0d0
      t_gbet = 0.0d0
      t_viol = 0.0d0
      t_gviol = 0.0d0
      n_func = 0
      n_grad = 0
      n_fhel = 0
      n_fbet = 0
      n_ghel = 0
      n_gbet = 0
      n_viol = 0
      n_gviol = 0
      n_map = 0
#ifndef SPLITELE
      nprint_ene=nprint_ene-1
#endif
      return
      end
c-------------------------------------------------------------------------
      block data nazwy
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'sizesclu.dat'
      include 'COMMON.NAMES'
      include 'COMMON.FFIELD'
      data restyp /
     &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
     & 'DSG','DGN','DSN','DTH',
     &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
     &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
     &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
     &'AIB','ABU','D'/
      data onelet /
     &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
     &'a','y','w','v','l','i','f','m','c','x',
     &'C','M','F','I','L','V','W','Y','A','G','T',
     &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
      data potname /'LJ','LJK','BP','GB','GBV'/
      data ename /
     &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
     &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
     &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
     &   "ESTR","ESCCOR","EVDW2_14","","EVDW_T"/
      data wname /
     &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
     &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
     &   "WHPB","WVDWPP","WBOND","WSCCOR","WSCP14","","WSC"/
      data nprint_ene /21/
      data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
     &  16,15,17,20,21/
      end 
c---------------------------------------------------------------------------
      subroutine init_int_table
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'sizesclu.dat'
#ifdef MPI
      include 'mpif.h'
#endif
#ifdef MPL
      include 'COMMON.INFO'
#endif
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.LOCAL'
      include 'COMMON.SBRIDGE'
      include 'COMMON.IOUNITS'
      logical scheck,lprint
#ifdef MPL
      integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
     & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
C... Determine the numbers of start and end SC-SC interaction 
C... to deal with by current processor.
      lprint=.false.
      if (lprint)
     &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
      n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
      MyRank=MyID-(MyGroup-1)*fgProcs
      call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
      if (lprint)
     &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
     &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
     &  ' my_sc_inde',my_sc_inde
      ind_sctint=0
      iatsc_s=0
      iatsc_e=0
#endif
      lprint=.false.
      do i=1,maxres
        nint_gr(i)=0
        nscp_gr(i)=0
        do j=1,maxint_gr
          istart(i,1)=0
          iend(i,1)=0
          ielstart(i)=0
          ielend(i)=0
          iscpstart(i,1)=0
          iscpend(i,1)=0    
        enddo
      enddo
      ind_scint=0
      ind_scint_old=0
cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
cd   &   (ihpb(i),jhpb(i),i=1,nss)
      do i=nnt,nct-1
        scheck=.false.
        do ii=1,nss
          if (ihpb(ii).eq.i+nres) then
            scheck=.true.
            jj=jhpb(ii)-nres
            goto 10
          endif
        enddo
   10   continue
cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
        if (scheck) then
          if (jj.eq.i+1) then
#ifdef MPL
            write (iout,*) 'jj=i+1'
            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
     & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
#else
            nint_gr(i)=1
            istart(i,1)=i+2
            iend(i,1)=nct
#endif
          else if (jj.eq.nct) then
#ifdef MPL
            write (iout,*) 'jj=nct'
            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
     &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
#else
            nint_gr(i)=1
            istart(i,1)=i+1
            iend(i,1)=nct-1
#endif
          else
#ifdef MPL
            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
     & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
            ii=nint_gr(i)+1
            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
     & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
#else
            nint_gr(i)=2
            istart(i,1)=i+1
            iend(i,1)=jj-1
            istart(i,2)=jj+1
            iend(i,2)=nct
#endif
          endif
        else
#ifdef MPL
          call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
     &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
#else
          nint_gr(i)=1
          istart(i,1)=i+1
          iend(i,1)=nct
          ind_scint=int_scint+nct-i
#endif
        endif
#ifdef MPL
        ind_scint_old=ind_scint
#endif
      enddo
   12 continue
#ifndef MPL
      iatsc_s=nnt
      iatsc_e=nct-1
#endif
#ifdef MPL
      if (lprint) then
        write (iout,*) 'Processor',MyID,' Group',MyGroup
        write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
      endif
#endif
      if (lprint) then
      write (iout,'(a)') 'Interaction array:'
      do i=iatsc_s,iatsc_e
        write (iout,'(i3,2(2x,2i3))') 
     & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
      enddo
      endif
      ispp=2
#ifdef MPL
C Now partition the electrostatic-interaction array
      npept=nct-nnt
      nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
      call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
      if (lprint)
     & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
     &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
     &               ' my_ele_inde',my_ele_inde
      iatel_s=0
      iatel_e=0
      ind_eleint=0
      ind_eleint_old=0
      do i=nnt,nct-3
        ijunk=0
        call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
     &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
      enddo ! i 
   13 continue
#else
      iatel_s=nnt
      iatel_e=nct-3
      do i=iatel_s,iatel_e
        ielstart(i)=i+2
        ielend(i)=nct-1
      enddo
#endif
      if (lprint) then
        write (iout,'(a)') 'Electrostatic interaction array:'
        do i=iatel_s,iatel_e
          write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
        enddo
      endif ! lprint
c     iscp=3
      iscp=2
C Partition the SC-p interaction array
#ifdef MPL
      nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
      call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
      if (lprint)
     & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
     &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
     &               ' my_scp_inde',my_scp_inde
      iatscp_s=0
      iatscp_e=0
      ind_scpint=0
      ind_scpint_old=0
      do i=nnt,nct-1
        if (i.lt.nnt+iscp) then
cd        write (iout,*) 'i.le.nnt+iscp'
          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
     &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
     &      iscpend(i,1),*14)
        else if (i.gt.nct-iscp) then
cd        write (iout,*) 'i.gt.nct-iscp'
          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
     &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
     &      iscpend(i,1),*14)
        else
          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
     &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
     &      iscpend(i,1),*14)
          ii=nscp_gr(i)+1
          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
     &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
     &      iscpend(i,ii),*14)
        endif
      enddo ! i
   14 continue
#else
      iatscp_s=nnt
      iatscp_e=nct-1
      do i=nnt,nct-1
        if (i.lt.nnt+iscp) then
          nscp_gr(i)=1
          iscpstart(i,1)=i+iscp
          iscpend(i,1)=nct
        elseif (i.gt.nct-iscp) then
          nscp_gr(i)=1
          iscpstart(i,1)=nnt
          iscpend(i,1)=i-iscp
        else
          nscp_gr(i)=2
          iscpstart(i,1)=nnt
          iscpend(i,1)=i-iscp
          iscpstart(i,2)=i+iscp
          iscpend(i,2)=nct
        endif 
      enddo ! i
#endif
      if (lprint) then
        write (iout,'(a)') 'SC-p interaction array:'
        do i=iatscp_s,iatscp_e
          write (iout,'(i3,2(2x,2i3))') 
     &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
        enddo
      endif ! lprint
C Partition local interactions
#ifdef MPL
      call int_bounds(nres-2,loc_start,loc_end)
      loc_start=loc_start+1
      loc_end=loc_end+1
      call int_bounds(nres-2,ithet_start,ithet_end)
      ithet_start=ithet_start+2
      ithet_end=ithet_end+2
      call int_bounds(nct-nnt-2,iphi_start,iphi_end) 
      iphi_start=iphi_start+nnt+2
      iphi_end=iphi_end+nnt+2
      call int_bounds(nres-3,itau_start,itau_end)
      itau_start=itau_start+3
      itau_end=itau_end+3
      if (lprint) then 
        write (iout,*) 'Processor:',MyID,
     & ' loc_start',loc_start,' loc_end',loc_end,
     & ' ithet_start',ithet_start,' ithet_end',ithet_end,
     & ' iphi_start',iphi_start,' iphi_end',iphi_end
        write (*,*) 'Processor:',MyID,
     & ' loc_start',loc_start,' loc_end',loc_end,
     & ' ithet_start',ithet_start,' ithet_end',ithet_end,
     & ' iphi_start',iphi_start,' iphi_end',iphi_end
      endif
      if (fgprocs.gt.1 .and. MyID.eq.BossID) then
        write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
     & nele_int_tot,' electrostatic and ',nscp_int_tot,
     & ' SC-p interactions','were distributed among',fgprocs,
     & ' fine-grain processors.'
      endif
#else
      loc_start=2
      loc_end=nres-1
      ithet_start=3 
      ithet_end=nres
      iphi_start=nnt+3
      iphi_end=nct
      itau_start=4
      itau_end=nres
#endif
      return
      end 
c---------------------------------------------------------------------------
      subroutine int_partition(int_index,lower_index,upper_index,atom,
     & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      integer int_index,lower_index,upper_index,atom,at_start,at_end,
     & first_atom,last_atom,int_gr,jat_start,jat_end
      logical lprn
      lprn=.false.
      if (lprn) write (iout,*) 'int_index=',int_index
      int_index_old=int_index
      int_index=int_index+last_atom-first_atom+1
      if (lprn) 
     &   write (iout,*) 'int_index=',int_index,
     &               ' int_index_old',int_index_old,
     &               ' lower_index=',lower_index,
     &               ' upper_index=',upper_index,
     &               ' atom=',atom,' first_atom=',first_atom,
     &               ' last_atom=',last_atom
      if (int_index.ge.lower_index) then
        int_gr=int_gr+1
        if (at_start.eq.0) then
          at_start=atom
          jat_start=first_atom-1+lower_index-int_index_old
        else
          jat_start=first_atom
        endif
        if (lprn) write (iout,*) 'jat_start',jat_start
        if (int_index.ge.upper_index) then
          at_end=atom
          jat_end=first_atom-1+upper_index-int_index_old
          return1
        else
          jat_end=last_atom
        endif
        if (lprn) write (iout,*) 'jat_end',jat_end
      endif
      return
      end
c------------------------------------------------------------------------------
      subroutine hpb_partition
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.SBRIDGE'
      include 'COMMON.IOUNITS'
#ifdef MPL
      include 'COMMON.INFO'
      call int_bounds(nhpb,link_start,link_end)
#else
      link_start=1
      link_end=nhpb
#endif
cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
cd   &  ' nhpb',nhpb,' link_start=',link_start,
cd   &  ' link_end',link_end
      return
      end