source/unres/src_MD/src/dihed_cons.F

   1       subroutine secstrp2dihc
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.GEO'
   5       include 'COMMON.BOUNDS'
   6       include 'COMMON.CHAIN'
   7       include 'COMMON.TORCNSTR'
   8         include 'COMMON.IOUNITS'
   9       character*1 secstruc(maxres)
  10       COMMON/SECONDARYS/secstruc
  11       character*80 line
  12       logical errflag
  13       external ilen
  14
  15 cdr      call getenv_loc('SECPREDFIL',secpred)
  16       lenpre=ilen(prefix)
  17       secpred=prefix(:lenpre)//'.spred'
  18
  19 #if defined(WINIFL) || defined(WINPGI)
  20       open(isecpred,file=secpred,status='old',readonly,shared)
  21 #elif (defined CRAY) || (defined AIX)
  22       open(isecpred,file=secpred,status='old',action='read')
  23 #elif (defined G77)
  24       open(isecpred,file=secpred,status='old')
  25 #else
  26       open(isecpred,file=secpred,status='old',readonly)
  27 #endif
  28 C read secondary structure prediction from JPRED here!
  29 !      read(isecpred,'(A80)',err=100,end=100) line
  30 !      read(line,'(f10.3)',err=110) ftors
  31        read(isecpred,'(f10.3)',err=110) ftors
  32
  33       write (iout,*) 'FTORS factor =',ftors
  34 ! initialize secstruc to any
  35        do i=1,nres
  36         secstruc(i) ='-'
  37       enddo
  38       ndih_constr=0
  39       ndih_nconstr=0
  40
  41       call read_secstr_pred(isecpred,iout,errflag)
  42       if (errflag) then
  43          write(iout,*)'There is a problem with the list of secondary-',
  44      &     'structure prediction'
  45          goto 100
  46       endif
  47 C 8/13/98 Set limits to generating the dihedral angles
  48       do i=1,nres
  49         phibound(1,i)=-pi
  50         phibound(2,i)=pi
  51       enddo
  52
  53       ii=0
  54       do i=1,nres
  55         if ( secstruc(i) .eq. 'H') then
  56 C Helix restraints for this residue
  57            ii=ii+1
  58            idih_constr(ii)=i
  59            phi0(ii) = 45.0D0*deg2rad
  60            drange(ii)= 5.0D0*deg2rad
  61            phibound(1,i) = phi0(ii)-drange(ii)
  62            phibound(2,i) = phi0(ii)+drange(ii)
  63         else if (secstruc(i) .eq. 'E') then
  64 C strand restraints for this residue
  65            ii=ii+1
  66            idih_constr(ii)=i
  67            phi0(ii) = 180.0D0*deg2rad
  68            drange(ii)= 5.0D0*deg2rad
  69            phibound(1,i) = phi0(ii)-drange(ii)
  70            phibound(2,i) = phi0(ii)+drange(ii)
  71         else
  72 C no restraints for this residue
  73            ndih_nconstr=ndih_nconstr+1
  74            idih_nconstr(ndih_nconstr)=i
  75         endif
  76       enddo
  77       ndih_constr=ii
  78       return
  79 100   continue
  80       write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
  81       return
  82  110  continue
  83       write(iout,'(A20)')'Error reading FTORS'
  84       return
  85       end
  86
  87       subroutine read_secstr_pred(jin,jout,errors)
  88
  89       implicit real*8 (a-h,o-z)
  90       INCLUDE 'DIMENSIONS'
  91       include 'COMMON.IOUNITS'
  92       include 'COMMON.CHAIN'
  93       character*1 secstruc(maxres)
  94       COMMON/SECONDARYS/secstruc
  95       EXTERNAL ILEN
  96       character*80 line,line1,ucase
  97       logical errflag,errors,blankline
  98
  99       errors=.false.
 100       read (jin,'(a)') line
 101       write (jout,'(2a)') '> ',line(1:78)
 102       line1=ucase(line)
 103 C Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
 104 C correspond to the end-groups.  ADD to the secondary structure prediction "-" for the
 105 C end-groups in the input file "*.spred"
 106
 107       iseq=1
 108       do while (index(line1,'$END').eq.0)
 109 * Override commented lines.
 110          ipos=1
 111          blankline=.false.
 112          do while (.not.blankline)
 113             line1=' '
 114             call mykey(line,line1,ipos,blankline,errflag)
 115             if (errflag) write (jout,'(2a)')
 116      & 'Error when reading sequence in line: ',line
 117             errors=errors .or. errflag
 118             if (.not. blankline .and. .not. errflag) then
 119                ipos1=2
 120                iend=ilen(line1)
 121                if (iseq.le.maxres) then
 122                   if (line1(1:1).eq.'-' ) then
 123                      secstruc(iseq)=line1(1:1)
 124                   else if ( ( ucase(line1(1:1)).eq.'E' ) .or.
 125      &                      ( ucase(line1(1:1)).eq.'H' ) ) then
 126                      secstruc(iseq)=ucase(line1(1:1))
 127                   else
 128                      errors=.true.
 129                      write (jout,1010) line1(1:1), iseq
 130                      goto 80
 131                   endif
 132                else
 133                   errors=.true.
 134                   write (jout,1000) iseq,maxres
 135                   goto 80
 136                endif
 137                do while (ipos1.le.iend)
 138
 139                   iseq=iseq+1
 140                   il=1
 141                   ipos1=ipos1+1
 142                   if (iseq.le.maxres) then
 143                      if (line1(ipos1-1:ipos1-1).eq.'-' ) then
 144                         secstruc(iseq)=line1(ipos1-1:ipos1-1)
 145                      else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or.
 146      &                     (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
 147                         secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
 148                      else
 149                         errors=.true.
 150                         write (jout,1010) line1(ipos1-1:ipos1-1), iseq
 151                         goto 80
 152                      endif
 153                   else
 154                      errors=.true.
 155                      write (jout,1000) iseq,maxres
 156                      goto 80
 157                   endif
 158                enddo
 159                iseq=iseq+1
 160             endif
 161          enddo
 162          read (jin,'(a)') line
 163          write (jout,'(2a)') '> ',line(1:78)
 164          line1=ucase(line)
 165       enddo
 166
 167 cd    write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
 168
 169 cd check whether the found length of the chain is correct.
 170       length_of_chain=iseq-1
 171       if (length_of_chain .ne. nres) then
 172 !        errors=.true.
 173         write (jout,'(a,i4,a,i4,a)')
 174      & 'Error: the number of labels specified in $SEC_STRUC_PRED ('
 175      & ,length_of_chain,') does not match with the number of residues ('
 176      & ,nres,').'
 177       endif
 178    80 continue
 179
 180  1000 format('Error - the number of residues (',i4,
 181      & ') has exceeded maximum (',i4,').')
 182  1010 format ('Error - unrecognized secondary structure label',a4,
 183      & ' in position',i4)
 184       return
 185       end