Merge branch 'bartek' of mmka.chem.univ.gda.pl:unres into bartek
[unres.git] / source / unres / src_MD / 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',action='read')
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