Fixed the homology gradient in finegrain mode
[unres.git] / source / wham / src-M / molread_zs.F
1       subroutine molread(*)
2 C
3 C Read molecular data.
4 C
5       implicit real*8 (a-h,o-z)
6       include 'DIMENSIONS'
7       include 'DIMENSIONS.ZSCOPT'
8       include 'COMMON.IOUNITS'
9       include 'COMMON.GEO'
10       include 'COMMON.VAR'
11       include 'COMMON.INTERACT'
12       include 'COMMON.LOCAL'
13       include 'COMMON.NAMES'
14       include 'COMMON.CHAIN'
15       include 'COMMON.FFIELD'
16       include 'COMMON.SBRIDGE'
17       include 'COMMON.TORCNSTR'
18       include 'COMMON.CONTROL'
19       character*4 sequence(maxres)
20       integer rescode
21       double precision x(maxvar)
22       character*320 controlcard,ucase
23       dimension itype_pdb(maxres)
24       logical seq_comp
25       call card_concat(controlcard,.true.)
26       call reada(controlcard,'SCAL14',scal14,0.4d0)
27       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
28       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
29       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
30       r0_corr=cutoff_corr-delt_corr
31       call readi(controlcard,"NRES",nres,0)
32       iscode=index(controlcard,"ONE_LETTER")
33       if (nres.le.0) then
34         write (iout,*) "Error: no residues in molecule"
35         return1
36       endif
37       if (nres.gt.maxres) then
38         write (iout,*) "Error: too many residues",nres,maxres
39       endif
40       write(iout,*) 'nres=',nres
41 C Read sequence of the protein
42       if (iscode.gt.0) then
43         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
44       else
45         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
46       endif
47 C Convert sequence to numeric code
48       do i=1,nres
49         itype(i)=rescode(i,sequence(i),iscode)
50       enddo
51       write (iout,*) "Numeric code:"
52       write (iout,'(20i4)') (itype(i),i=1,nres)
53       do i=1,nres-1
54 #ifdef PROCOR
55         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
56 #else
57         if (itype(i).eq.21) then
58 #endif
59           itel(i)=0
60 #ifdef PROCOR
61         else if (itype(i+1).ne.20) then
62 #else
63         else if (itype(i).ne.20) then
64 #endif
65           itel(i)=1
66         else
67           itel(i)=2
68         endif
69       enddo
70       call read_bridge
71
72       if (with_dihed_constr) then
73
74       read (inp,*) ndih_constr
75       if (ndih_constr.gt.0) then
76         read (inp,*) ftors
77         write (iout,*) 'FTORS',ftors
78         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
79         write (iout,*)
80      &   'There are',ndih_constr,' constraints on phi angles.'
81         do i=1,ndih_constr
82           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
83         enddo
84         do i=1,ndih_constr
85           phi0(i)=deg2rad*phi0(i)
86           drange(i)=deg2rad*drange(i)
87         enddo
88       endif
89
90       endif
91
92       nnt=1
93       nct=nres
94       if (itype(1).eq.21) nnt=2
95       if (itype(nres).eq.21) nct=nct-1
96       write(iout,*) 'NNT=',NNT,' NCT=',NCT
97       call setup_var
98       call init_int_table
99       if (ns.gt.0) then
100         write (iout,'(/a,i3,a)') 'The chain contains',ns,
101      &  ' disulfide-bridging cysteines.'
102         write (iout,'(20i4)') (iss(i),i=1,ns)
103         write (iout,'(/a/)') 'Pre-formed links are:' 
104         do i=1,nss
105           i1=ihpb(i)-nres
106           i2=jhpb(i)-nres
107           it1=itype(i1)
108           it2=itype(i2)
109          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
110      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
111      &    dhpb(i),ebr,forcon(i)
112         enddo
113       endif
114       write (iout,'(a)')
115       return
116       end
117 c-----------------------------------------------------------------------------
118       logical function seq_comp(itypea,itypeb,length)
119       implicit none
120       integer length,itypea(length),itypeb(length)
121       integer i
122       do i=1,length
123         if (itypea(i).ne.itypeb(i)) then
124           seq_comp=.false.
125           return
126         endif
127       enddo
128       seq_comp=.true.
129       return
130       end
131 c-----------------------------------------------------------------------------
132       subroutine read_bridge
133 C Read information about disulfide bridges.
134       implicit real*8 (a-h,o-z)
135       include 'DIMENSIONS'
136       include 'DIMENSIONS.ZSCOPT'
137       include 'COMMON.IOUNITS'
138       include 'COMMON.GEO'
139       include 'COMMON.VAR'
140       include 'COMMON.INTERACT'
141       include 'COMMON.NAMES'
142       include 'COMMON.CHAIN'
143       include 'COMMON.FFIELD'
144       include 'COMMON.SBRIDGE'
145 C Read bridging residues.
146       read (inp,*) ns,(iss(i),i=1,ns)
147       print *,'ns=',ns
148       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
149 C Check whether the specified bridging residues are cystines.
150       do i=1,ns
151         if (itype(iss(i)).ne.1) then
152           write (iout,'(2a,i3,a)') 
153      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
154      &   ' can form a disulfide bridge?!!!'
155           write (*,'(2a,i3,a)') 
156      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
157      &   ' can form a disulfide bridge?!!!'
158          stop
159         endif
160       enddo
161 C Read preformed bridges.
162       if (ns.gt.0) then
163       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
164       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
165       if (nss.gt.0) then
166         nhpb=nss
167 C Check if the residues involved in bridges are in the specified list of
168 C bridging residues.
169         do i=1,nss
170           do j=1,i-1
171             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
172      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
173               write (iout,'(a,i3,a)') 'Disulfide pair',i,
174      &      ' contains residues present in other pairs.'
175               write (*,'(a,i3,a)') 'Disulfide pair',i,
176      &      ' contains residues present in other pairs.'
177               stop 
178             endif
179           enddo
180           do j=1,ns
181             if (ihpb(i).eq.iss(j)) goto 10
182           enddo
183           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
184    10     continue
185           do j=1,ns
186             if (jhpb(i).eq.iss(j)) goto 20
187           enddo
188           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
189    20     continue
190           dhpb(i)=dbr
191           forcon(i)=fbr
192         enddo
193         do i=1,nss
194           ihpb(i)=ihpb(i)+nres
195           jhpb(i)=jhpb(i)+nres
196         enddo
197       endif
198       endif
199       return
200       end
201 c------------------------------------------------------------------------------
202       subroutine read_angles(kanal,iscor,energ,iprot,*)
203       implicit real*8 (a-h,o-z)
204       include 'DIMENSIONS'
205       include 'DIMENSIONS.ZSCOPT'
206       include 'COMMON.INTERACT'
207       include 'COMMON.SBRIDGE'
208       include 'COMMON.GEO'
209       include 'COMMON.VAR'
210       include 'COMMON.CHAIN'
211       include 'COMMON.IOUNITS'
212       character*80 lineh
213       read(kanal,'(a80)',end=10,err=10) lineh
214       read(lineh(:5),*,err=8) ic
215       read(lineh(6:),*,err=8) energ
216       goto 9
217     8 ic=1
218       print *,'error, assuming e=1d10',lineh
219       energ=1d10
220       nss=0
221     9 continue
222       read(lineh(18:),*,end=10,err=10) nss
223       IF (NSS.LT.9) THEN
224         read (lineh(20:),*,end=10,err=10)
225      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
226       ELSE
227         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
228         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
229      &    I=9,NSS),iscor
230       ENDIF
231 c      print *,"energy",energ," iscor",iscor
232       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
233       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
234       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
235       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
236       do i=1,nres
237         theta(i)=deg2rad*theta(i)
238         phi(i)=deg2rad*phi(i)
239         alph(i)=deg2rad*alph(i)
240         omeg(i)=deg2rad*omeg(i)
241       enddo
242       return
243    10 return1
244       end