changes in wham
[unres.git] / source / unres / src_MD-M / 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(1)
32
33       write (iout,*) 'FTORS factor =',ftors(1)
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          ftors(i)=ftors(1)
56         if ( secstruc(i) .eq. 'H') then
57 C Helix restraints for this residue               
58            ii=ii+1 
59            idih_constr(ii)=i
60            phi0(ii) = 45.0D0*deg2rad
61            drange(ii)= 5.0D0*deg2rad
62            phibound(1,i) = phi0(ii)-drange(ii)
63            phibound(2,i) = phi0(ii)+drange(ii)
64         else if (secstruc(i) .eq. 'E') then
65 C strand restraints for this residue               
66            ii=ii+1 
67            idih_constr(ii)=i 
68            phi0(ii) = 180.0D0*deg2rad
69            drange(ii)= 5.0D0*deg2rad
70            phibound(1,i) = phi0(ii)-drange(ii)
71            phibound(2,i) = phi0(ii)+drange(ii)
72         else
73 C no restraints for this residue               
74            ndih_nconstr=ndih_nconstr+1
75            idih_nconstr(ndih_nconstr)=i
76         endif        
77       enddo
78       ndih_constr=ii
79       return
80 100   continue
81       write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
82       return
83  110  continue
84       write(iout,'(A20)')'Error reading FTORS'
85       return
86       end 
87
88       subroutine read_secstr_pred(jin,jout,errors)
89
90       implicit real*8 (a-h,o-z)
91       INCLUDE 'DIMENSIONS'
92       include 'COMMON.IOUNITS'
93       include 'COMMON.CHAIN'
94       character*1 secstruc(maxres)
95       COMMON/SECONDARYS/secstruc
96       EXTERNAL ILEN
97       character*80 line,line1,ucase
98       logical errflag,errors,blankline
99
100       errors=.false.
101       read (jin,'(a)') line
102       write (jout,'(2a)') '> ',line(1:78)
103       line1=ucase(line)
104 C Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
105 C correspond to the end-groups.  ADD to the secondary structure prediction "-" for the
106 C end-groups in the input file "*.spred"
107
108       iseq=1
109       do while (index(line1,'$END').eq.0)
110 * Override commented lines.
111          ipos=1
112          blankline=.false.
113          do while (.not.blankline)
114             line1=' '
115             call mykey(line,line1,ipos,blankline,errflag) 
116             if (errflag) write (jout,'(2a)')
117      & 'Error when reading sequence in line: ',line
118             errors=errors .or. errflag
119             if (.not. blankline .and. .not. errflag) then
120                ipos1=2
121                iend=ilen(line1)
122                if (iseq.le.maxres) then
123                   if (line1(1:1).eq.'-' ) then 
124                      secstruc(iseq)=line1(1:1)
125                   else if ( ( ucase(line1(1:1)).eq.'E' ) .or.
126      &                      ( ucase(line1(1:1)).eq.'H' ) ) then
127                      secstruc(iseq)=ucase(line1(1:1))
128                   else
129                      errors=.true.
130                      write (jout,1010) line1(1:1), iseq
131                      goto 80
132                   endif                  
133                else
134                   errors=.true.
135                   write (jout,1000) iseq,maxres
136                   goto 80
137                endif
138                do while (ipos1.le.iend)
139
140                   iseq=iseq+1
141                   il=1
142                   ipos1=ipos1+1
143                   if (iseq.le.maxres) then
144                      if (line1(ipos1-1:ipos1-1).eq.'-' ) then 
145                         secstruc(iseq)=line1(ipos1-1:ipos1-1)
146                      else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or.
147      &                     (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
148                         secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
149                      else
150                         errors=.true.
151                         write (jout,1010) line1(ipos1-1:ipos1-1), iseq
152                         goto 80
153                      endif                  
154                   else
155                      errors=.true.
156                      write (jout,1000) iseq,maxres
157                      goto 80
158                   endif
159                enddo
160                iseq=iseq+1
161             endif
162          enddo
163          read (jin,'(a)') line
164          write (jout,'(2a)') '> ',line(1:78)
165          line1=ucase(line)
166       enddo
167
168 cd    write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
169
170 cd check whether the found length of the chain is correct.
171       length_of_chain=iseq-1
172       if (length_of_chain .ne. nres) then
173 !        errors=.true.
174         write (jout,'(a,i4,a,i4,a)') 
175      & 'Error: the number of labels specified in $SEC_STRUC_PRED ('
176      & ,length_of_chain,') does not match with the number of residues ('
177      & ,nres,').' 
178       endif
179    80 continue
180
181  1000 format('Error - the number of residues (',i4,
182      & ') has exceeded maximum (',i4,').')
183  1010 format ('Error - unrecognized secondary structure label',a4,
184      & ' in position',i4)
185       return
186       end