subroutine molread(*)
C
C Read molecular data.
C
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.LOCAL'
      include 'COMMON.NAMES'
      include 'COMMON.CHAIN'
      include 'COMMON.FFIELD'
      include 'COMMON.SBRIDGE'
      include 'COMMON.TORCNSTR'
      include 'COMMON.CONTROL'
      character*4 sequence(maxres)
      integer rescode
      double precision x(maxvar)
      character*320 controlcard,ucase
      dimension itype_pdb(maxres)
      logical seq_comp
      call card_concat(controlcard,.true.)
      call reada(controlcard,'SCAL14',scal14,0.4d0)
      call reada(controlcard,'SCALSCP',scalscp,1.0d0)
      call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
      call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
C     Bartek
      call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
      call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
      call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
      call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
      write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
     &  " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
      r0_corr=cutoff_corr-delt_corr
      call readi(controlcard,"NRES",nres,0)
      iscode=index(controlcard,"ONE_LETTER")
      if (nres.le.0) then
        write (iout,*) "Error: no residues in molecule"
        return1
      endif
      if (nres.gt.maxres) then
        write (iout,*) "Error: too many residues",nres,maxres
      endif
      write(iout,*) 'nres=',nres
C Read sequence of the protein
      if (iscode.gt.0) then
        read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
      else
        read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
      endif
C Convert sequence to numeric code
      do i=1,nres
        itype(i)=rescode(i,sequence(i),iscode)
      enddo
      write (iout,*) "Numeric code:"
      write (iout,'(20i4)') (itype(i),i=1,nres)
      do i=1,nres-1
#ifdef PROCOR
        if (itype(i).eq.21 .or. itype(i+1).eq.21) then
#else
        if (itype(i).eq.21) then
#endif
          itel(i)=0
#ifdef PROCOR
        else if (itype(i+1).ne.20) then
#else
        else if (itype(i).ne.20) then
#endif
          itel(i)=1
        else
          itel(i)=2
        endif
      enddo
      call read_bridge

      if (with_dihed_constr) then

      read (inp,*) ndih_constr
      if (ndih_constr.gt.0) then
        read (inp,*) ftors
        write (iout,*) 'FTORS',ftors
        read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
        write (iout,*)
     &   'There are',ndih_constr,' constraints on phi angles.'
        do i=1,ndih_constr
          write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
        enddo
        do i=1,ndih_constr
          phi0(i)=deg2rad*phi0(i)
          drange(i)=deg2rad*drange(i)
        enddo
      endif

      endif

      nnt=1
      nct=nres
      if (itype(1).eq.21) nnt=2
      if (itype(nres).eq.21) nct=nct-1
      write(iout,*) 'NNT=',NNT,' NCT=',NCT

C     Juyong:READ init_vars
C     Initialize variables!
C     Juyong:READ read_info
C     READ fragment information!!
C     both routines should be in dfa.F file!!

      if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
     &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
       write (iout,*) "Calling init_dfa_vars"
       call flush(iout)
       call init_dfa_vars
       write (iout,*) 'init_dfa_vars finished!'
       call flush(iout)
       call read_dfa_info
       write (iout,*) 'read_dfa_info finished!'
       call flush(iout)
      endif

c Read distance restraints
      if (constr_dist.gt.0) then
        if (refstr) call read_ref_structure(*11)
        call read_dist_constr
        call hpb_partition
      endif

      if (constr_homology.gt.0) then
        call read_constr_homology
      endif


      call setup_var
      call init_int_table
      if (ns.gt.0) then
        write (iout,'(/a,i3,a)') 'The chain contains',ns,
     &  ' disulfide-bridging cysteines.'
        write (iout,'(20i4)') (iss(i),i=1,ns)
        write (iout,'(/a/)') 'Pre-formed links are:' 
        do i=1,nss
	  i1=ihpb(i)-nres
	  i2=jhpb(i)-nres
	  it1=itype(i1)
	  it2=itype(i2)
         write (iout,'(2a,i3,3a,i3,a,3f10.3)')
     &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
     &    dhpb(i),ebr,forcon(i)
        enddo
      endif
      write (iout,'(a)')
      return
   11 stop "Error reading reference structure"
      end
c-----------------------------------------------------------------------------
      logical function seq_comp(itypea,itypeb,length)
      implicit none
      integer length,itypea(length),itypeb(length)
      integer i
      do i=1,length
        if (itypea(i).ne.itypeb(i)) then
          seq_comp=.false.
          return
        endif
      enddo
      seq_comp=.true.
      return
      end
c-----------------------------------------------------------------------------
      subroutine read_bridge
C Read information about disulfide bridges.
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.NAMES'
      include 'COMMON.CHAIN'
      include 'COMMON.FFIELD'
      include 'COMMON.SBRIDGE'
C Read bridging residues.
      read (inp,*) ns,(iss(i),i=1,ns)
      print *,'ns=',ns
      write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
C Check whether the specified bridging residues are cystines.
      do i=1,ns
	if (itype(iss(i)).ne.1) then
	  write (iout,'(2a,i3,a)') 
     &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
     &   ' can form a disulfide bridge?!!!'
	  write (*,'(2a,i3,a)') 
     &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
     &   ' can form a disulfide bridge?!!!'
         stop
        endif
      enddo
C Read preformed bridges.
      if (ns.gt.0) then
      read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
      write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
      if (nss.gt.0) then
        nhpb=nss
C Check if the residues involved in bridges are in the specified list of
C bridging residues.
        do i=1,nss
          do j=1,i-1
	    if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
     &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
	      write (iout,'(a,i3,a)') 'Disulfide pair',i,
     &      ' contains residues present in other pairs.'
	      write (*,'(a,i3,a)') 'Disulfide pair',i,
     &      ' contains residues present in other pairs.'
              stop 
	    endif
          enddo
	  do j=1,ns
	    if (ihpb(i).eq.iss(j)) goto 10
          enddo
          write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
   10     continue
	  do j=1,ns
	    if (jhpb(i).eq.iss(j)) goto 20
          enddo
          write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
   20     continue
          dhpb(i)=dbr
          forcon(i)=fbr
        enddo
        do i=1,nss
          ihpb(i)=ihpb(i)+nres
          jhpb(i)=jhpb(i)+nres
        enddo
      endif
      endif
      if (ns.gt.0.and.dyn_ss) then
C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
C the bond
          do i=nss+1,nhpb
C /06/28/2013 Adasko: nss number of full SS bonds
            ihpb(i-nss)=ihpb(i)
            jhpb(i-nss)=jhpb(i)
            forcon(i-nss)=forcon(i)
            dhpb(i-nss)=dhpb(i)
          enddo
          nhpb=nhpb-nss
          nss=0
          call hpb_partition
          do i=1,ns
            dyn_ss_mask(iss(i))=.true.
C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
          enddo
      endif
      return
      end
c------------------------------------------------------------------------------
      subroutine read_angles(kanal,iscor,energ,iprot,*)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'COMMON.INTERACT'
      include 'COMMON.SBRIDGE'
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      character*80 lineh
      read(kanal,'(a80)',end=10,err=10) lineh
      read(lineh(:5),*,err=8) ic
      read(lineh(6:),*,err=8) energ
      goto 9
    8 ic=1
      print *,'error, assuming e=1d10',lineh
      energ=1d10
      nss=0
    9 continue
      read(lineh(18:),*,end=10,err=10) nss
      IF (NSS.LT.9) THEN
        read (lineh(20:),*,end=10,err=10)
     &  (IHPB(I),JHPB(I),I=1,NSS),iscor
      ELSE
        read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
        read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
     &    I=9,NSS),iscor
      ENDIF
c      print *,"energy",energ," iscor",iscor
      read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
      read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
      read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
      read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
      do i=1,nres
        theta(i)=deg2rad*theta(i)
        phi(i)=deg2rad*phi(i)
        alph(i)=deg2rad*alph(i)
        omeg(i)=deg2rad*omeg(i)
      enddo
      return
   10 return1
      end
c-------------------------------------------------------------------------------
      subroutine read_dist_constr
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CONTROL'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      include 'COMMON.SBRIDGE'
      integer ifrag_(2,100),ipair_(2,100)
      double precision wfrag_(100),wpair_(100)
      character*500 controlcard
c      write (iout,*) "Calling read_dist_constr"
c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
c      call flush(iout)
      call card_concat(controlcard,.true.)
      call readi(controlcard,"NFRAG",nfrag_,0)
      call readi(controlcard,"NPAIR",npair_,0)
      call readi(controlcard,"NDIST",ndist_,0)
      call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
      call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
      call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
      call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
      call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
      write (iout,*) "IFRAG"
      do i=1,nfrag_
        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
      enddo
      write (iout,*) "IPAIR"
      do i=1,npair_
        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
      enddo
      call flush(iout)
      if (.not.refstr .and. nfrag_.gt.0) then
        write (iout,*) 
     &  "ERROR: no reference structure to compute distance restraints"
        write (iout,*)
     &  "Restraints must be specified explicitly (NDIST=number)"
        stop 
      endif
      if (nfrag_.lt.2 .and. npair_.gt.0) then 
        write (iout,*) "ERROR: Less than 2 fragments specified",
     &   " but distance restraints between pairs requested"
        stop 
      endif 
      call flush(iout)
      do i=1,nfrag_
        if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
        if (ifrag_(2,i).gt.nstart_sup+nsup-1)
     &    ifrag_(2,i)=nstart_sup+nsup-1
c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
        call flush(iout)
        if (wfrag_(i).gt.0.0d0) then
        do j=ifrag_(1,i),ifrag_(2,i)-1
          do k=j+1,ifrag_(2,i)
            write (iout,*) "j",j," k",k
            ddjk=dist(j,k)
            if (constr_dist.eq.1) then
            nhpb=nhpb+1
            ihpb(nhpb)=j
            jhpb(nhpb)=k
              dhpb(nhpb)=ddjk
            forcon(nhpb)=wfrag_(i) 
            else if (constr_dist.eq.2) then
              if (ddjk.le.dist_cut) then
                nhpb=nhpb+1
                ihpb(nhpb)=j
                jhpb(nhpb)=k
                dhpb(nhpb)=ddjk
                forcon(nhpb)=wfrag_(i) 
              endif
            else
              nhpb=nhpb+1
              ihpb(nhpb)=j
              jhpb(nhpb)=k
              dhpb(nhpb)=ddjk
              forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
            endif
            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
          enddo
        enddo
        endif
      enddo
      do i=1,npair_
        if (wpair_(i).gt.0.0d0) then
        ii = ipair_(1,i)
        jj = ipair_(2,i)
        if (ii.gt.jj) then
          itemp=ii
          ii=jj
          jj=itemp
        endif
        do j=ifrag_(1,ii),ifrag_(2,ii)
          do k=ifrag_(1,jj),ifrag_(2,jj)
            nhpb=nhpb+1
            ihpb(nhpb)=j
            jhpb(nhpb)=k
            forcon(nhpb)=wpair_(i)
            dhpb(nhpb)=dist(j,k)
            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
          enddo
        enddo
        endif
      enddo 
      do i=1,ndist_
        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
     &     ibecarb(i),forcon(nhpb+1)
        if (forcon(nhpb+1).gt.0.0d0) then
          nhpb=nhpb+1
          if (ibecarb(i).gt.0) then
            ihpb(i)=ihpb(i)+nres
            jhpb(i)=jhpb(i)+nres
          endif
          if (dhpb(nhpb).eq.0.0d0) 
     &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
        endif
      enddo
      do i=1,nhpb
          write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
     &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
      enddo
      call flush(iout)
      return
      end



c====-------------------------------------------------------------------

      subroutine read_constr_homology

      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.CONTROL'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      include 'COMMON.INTERACT'
      include 'COMMON.GEO'
      double precision odl_temp,sigma_odl_temp
      common /przechowalnia/ odl_temp(maxres,maxres,max_template),
     &    sigma_odl_temp(maxres,maxres,max_template)
      character*2 kic2
      character*24 model_ki_dist, model_ki_angle
      character*500 controlcard
      integer ki, i, j, k, l
      logical lprn /.true./

      call card_concat(controlcard,.true.)
      call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
      call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
      write (iout,*) "wga_dist",waga_dist," waga_angle",waga_angle

      lim_odl=0
      lim_dih=0
      do i=1,nres
        do j=i+2,nres
          do ki=1,constr_homology
            sigma_odl_temp(i,j,ki)=0.0d0
            odl_temp(i,j,ki)=0.0d0
          enddo
        enddo
      enddo
      do i=1,nres-3
        do ki=1,constr_homology
          dih(ki,i)=0.0d0
          sigma_dih(ki,i)=0.0d0
        enddo
      enddo
      do ki=1,constr_homology
          write(kic2,'(i2)') ki
          if (ki.le.9) kic2="0"//kic2(2:2)

          model_ki_dist="model"//kic2//".dist"
          model_ki_angle="model"//kic2//".angle"
        open (ientin,file=model_ki_dist,status='old')
        do irec=1,maxdim !petla do czytania wiezow na odleglosc
          read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
     &       sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
          odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
          sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
     &     sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
        enddo
 1401 continue
        close (ientin)
        open (ientin,file=model_ki_angle,status='old')
        do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
          read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
     &      sigma_dih(ki,i+nnt-1)
          if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
          sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
        enddo
 1402 continue
        close (ientin)
      enddo
      ii=0
      write (iout,*) "nnt",nnt," nct",nct
      do i=nnt,nct-2
        do j=i+2,nct
          ki=1
c          write (iout,*) "i",i," j",j," constr_homology",constr_homology
          do while (ki.le.constr_homology .and.
     &        sigma_odl_temp(i,j,ki).le.0.0d0)
c            write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
            ki=ki+1
          enddo
c          write (iout,*) "ki",ki
          if (ki.gt.constr_homology) cycle
          ii=ii+1
          ires_homo(ii)=i
          jres_homo(ii)=j
          do ki=1,constr_homology
            odl(ki,ii)=odl_temp(i,j,ki)
            sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
          enddo      
        enddo
      enddo
      lim_odl=ii
      if (constr_homology.gt.0) call homology_partition
c Print restraints
      if (.not.lprn) return
      write (iout,*) "Distance restraints from templates"
      do ii=1,lim_odl
        write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
     &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
      enddo
      write (iout,*) "Dihedral angle restraints from templates"
      do i=nnt,lim_dih
        write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
     &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
      enddo
c      write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
c      write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)

      return
      end