Adam's changes
[unres.git] / source / wham / src / 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       if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
52         write (iout,*)
53      &   "Glycine is the first full residue, initial dummy deleted"
54         do i=1,nres
55           itype(i)=itype(i+1)
56         enddo
57         nres=nres-1
58       endif
59       if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
60         write (iout,*)
61      &   "Glycine is the last full residue, terminal dummy deleted"
62         nres=nres-1
63       endif
64       write (iout,*) "Numeric code:"
65       write (iout,'(20i4)') (itype(i),i=1,nres)
66       do i=1,nres-1
67 #ifdef PROCOR
68         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
69 #else
70         if (itype(i).eq.21) then
71 #endif
72           itel(i)=0
73 #ifdef PROCOR
74         else if (itype(i+1).ne.20) then
75 #else
76         else if (itype(i).ne.20) then
77 #endif
78           itel(i)=1
79         else
80           itel(i)=2
81         endif
82       enddo
83       call read_bridge
84
85       if (with_dihed_constr) then
86
87       read (inp,*) ndih_constr
88       if (ndih_constr.gt.0) then
89         read (inp,*) ftors
90         write (iout,*) 'FTORS',ftors
91         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
92         write (iout,*)
93      &   'There are',ndih_constr,' constraints on phi angles.'
94         do i=1,ndih_constr
95           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
96         enddo
97         do i=1,ndih_constr
98           phi0(i)=deg2rad*phi0(i)
99           drange(i)=deg2rad*drange(i)
100         enddo
101       endif
102
103       endif
104
105       nnt=1
106       nct=nres
107       if (itype(1).eq.21) nnt=2
108       if (itype(nres).eq.21) nct=nct-1
109       write(iout,*) 'NNT=',NNT,' NCT=',NCT
110 c Read distance restraints
111       if (constr_dist.gt.0) then
112         if (refstr) call read_ref_structure(*11)
113         call read_dist_constr
114         call hpb_partition
115       endif
116
117       call setup_var
118       call init_int_table
119       if (ns.gt.0) then
120         write (iout,'(/a,i3,a)') 'The chain contains',ns,
121      &  ' disulfide-bridging cysteines.'
122         write (iout,'(20i4)') (iss(i),i=1,ns)
123         write (iout,'(/a/)') 'Pre-formed links are:' 
124         do i=1,nss
125           i1=ihpb(i)-nres
126           i2=jhpb(i)-nres
127           it1=itype(i1)
128           it2=itype(i2)
129          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
130      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
131      &    dhpb(i),ebr,forcon(i)
132         enddo
133       endif
134       write (iout,'(a)')
135       return
136    11 stop "Error reading reference structure"
137       end
138 c-----------------------------------------------------------------------------
139       logical function seq_comp(itypea,itypeb,length)
140       implicit none
141       integer length,itypea(length),itypeb(length)
142       integer i
143       do i=1,length
144         if (itypea(i).ne.itypeb(i)) then
145           seq_comp=.false.
146           return
147         endif
148       enddo
149       seq_comp=.true.
150       return
151       end
152 c-----------------------------------------------------------------------------
153       subroutine read_bridge
154 C Read information about disulfide bridges.
155       implicit real*8 (a-h,o-z)
156       include 'DIMENSIONS'
157       include 'DIMENSIONS.ZSCOPT'
158       include 'COMMON.IOUNITS'
159       include 'COMMON.GEO'
160       include 'COMMON.VAR'
161       include 'COMMON.INTERACT'
162       include 'COMMON.NAMES'
163       include 'COMMON.CHAIN'
164       include 'COMMON.FFIELD'
165       include 'COMMON.SBRIDGE'
166 C Read bridging residues.
167       read (inp,*) ns,(iss(i),i=1,ns)
168       print *,'ns=',ns
169       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
170 C Check whether the specified bridging residues are cystines.
171       do i=1,ns
172         if (itype(iss(i)).ne.1) then
173           write (iout,'(2a,i3,a)') 
174      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
175      &   ' can form a disulfide bridge?!!!'
176           write (*,'(2a,i3,a)') 
177      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
178      &   ' can form a disulfide bridge?!!!'
179          stop
180         endif
181       enddo
182 C Read preformed bridges.
183       if (ns.gt.0) then
184       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
185       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
186       if (nss.gt.0) then
187         nhpb=nss
188 C Check if the residues involved in bridges are in the specified list of
189 C bridging residues.
190         do i=1,nss
191           do j=1,i-1
192             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
193      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
194               write (iout,'(a,i3,a)') 'Disulfide pair',i,
195      &      ' contains residues present in other pairs.'
196               write (*,'(a,i3,a)') 'Disulfide pair',i,
197      &      ' contains residues present in other pairs.'
198               stop 
199             endif
200           enddo
201           do j=1,ns
202             if (ihpb(i).eq.iss(j)) goto 10
203           enddo
204           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
205    10     continue
206           do j=1,ns
207             if (jhpb(i).eq.iss(j)) goto 20
208           enddo
209           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
210    20     continue
211           dhpb(i)=dbr
212           forcon(i)=fbr
213         enddo
214         do i=1,nss
215           ihpb(i)=ihpb(i)+nres
216           jhpb(i)=jhpb(i)+nres
217         enddo
218       endif
219       endif
220       if (ns.gt.0.and.dyn_ss) then
221 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
222 C the bond
223           do i=nss+1,nhpb
224 C /06/28/2013 Adasko: nss number of full SS bonds
225             ihpb(i-nss)=ihpb(i)
226             jhpb(i-nss)=jhpb(i)
227             forcon(i-nss)=forcon(i)
228             dhpb(i-nss)=dhpb(i)
229           enddo
230           nhpb=nhpb-nss
231           nss=0
232           call hpb_partition
233           do i=1,ns
234             dyn_ss_mask(iss(i))=.true.
235 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
236 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
237           enddo
238       endif
239       return
240       end
241 c------------------------------------------------------------------------------
242       subroutine read_angles(kanal,iscor,energ,iprot,*)
243       implicit real*8 (a-h,o-z)
244       include 'DIMENSIONS'
245       include 'DIMENSIONS.ZSCOPT'
246       include 'COMMON.INTERACT'
247       include 'COMMON.SBRIDGE'
248       include 'COMMON.GEO'
249       include 'COMMON.VAR'
250       include 'COMMON.CHAIN'
251       include 'COMMON.IOUNITS'
252       character*80 lineh
253       read(kanal,'(a80)',end=10,err=10) lineh
254       read(lineh(:5),*,err=8) ic
255       read(lineh(6:),*,err=8) energ
256       goto 9
257     8 ic=1
258       print *,'error, assuming e=1d10',lineh
259       energ=1d10
260       nss=0
261     9 continue
262       read(lineh(18:),*,end=10,err=10) nss
263       IF (NSS.LT.9) THEN
264         read (lineh(20:),*,end=10,err=10)
265      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
266       ELSE
267         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
268         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
269      &    I=9,NSS),iscor
270       ENDIF
271 c      print *,"energy",energ," iscor",iscor
272       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
273       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
274       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
275       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
276       do i=1,nres
277         theta(i)=deg2rad*theta(i)
278         phi(i)=deg2rad*phi(i)
279         alph(i)=deg2rad*alph(i)
280         omeg(i)=deg2rad*omeg(i)
281       enddo
282       return
283    10 return1
284       end
285 c-------------------------------------------------------------------------------
286       subroutine read_dist_constr
287       implicit real*8 (a-h,o-z)
288       include 'DIMENSIONS'
289       include 'COMMON.CONTROL'
290       include 'COMMON.CHAIN'
291       include 'COMMON.IOUNITS'
292       include 'COMMON.SBRIDGE'
293       integer ifrag_(2,100),ipair_(2,100)
294       double precision wfrag_(100),wpair_(100)
295       character*500 controlcard
296 c      write (iout,*) "Calling read_dist_constr"
297 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
298 c      call flush(iout)
299       call card_concat(controlcard,.true.)
300       call readi(controlcard,"NFRAG",nfrag_,0)
301       call readi(controlcard,"NPAIR",npair_,0)
302       call readi(controlcard,"NDIST",ndist_,0)
303       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
304       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
305       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
306       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
307       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
308       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
309       write (iout,*) "IFRAG"
310       do i=1,nfrag_
311         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
312       enddo
313       write (iout,*) "IPAIR"
314       do i=1,npair_
315         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
316       enddo
317       call flush(iout)
318       if (.not.refstr .and. nfrag_.gt.0) then
319         write (iout,*) 
320      &  "ERROR: no reference structure to compute distance restraints"
321         write (iout,*)
322      &  "Restraints must be specified explicitly (NDIST=number)"
323         stop 
324       endif
325       if (nfrag_.lt.2 .and. npair_.gt.0) then 
326         write (iout,*) "ERROR: Less than 2 fragments specified",
327      &   " but distance restraints between pairs requested"
328         stop 
329       endif 
330       call flush(iout)
331       do i=1,nfrag_
332         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
333         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
334      &    ifrag_(2,i)=nstart_sup+nsup-1
335 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
336         call flush(iout)
337         if (wfrag_(i).gt.0.0d0) then
338         do j=ifrag_(1,i),ifrag_(2,i)-1
339           do k=j+1,ifrag_(2,i)
340             write (iout,*) "j",j," k",k
341             ddjk=dist(j,k)
342             if (constr_dist.eq.1) then
343             nhpb=nhpb+1
344             ihpb(nhpb)=j
345             jhpb(nhpb)=k
346               dhpb(nhpb)=ddjk
347             forcon(nhpb)=wfrag_(i) 
348             else if (constr_dist.eq.2) then
349               if (ddjk.le.dist_cut) then
350                 nhpb=nhpb+1
351                 ihpb(nhpb)=j
352                 jhpb(nhpb)=k
353                 dhpb(nhpb)=ddjk
354                 forcon(nhpb)=wfrag_(i) 
355               endif
356             else
357               nhpb=nhpb+1
358               ihpb(nhpb)=j
359               jhpb(nhpb)=k
360               dhpb(nhpb)=ddjk
361               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
362             endif
363             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
364      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
365           enddo
366         enddo
367         endif
368       enddo
369       do i=1,npair_
370         if (wpair_(i).gt.0.0d0) then
371         ii = ipair_(1,i)
372         jj = ipair_(2,i)
373         if (ii.gt.jj) then
374           itemp=ii
375           ii=jj
376           jj=itemp
377         endif
378         do j=ifrag_(1,ii),ifrag_(2,ii)
379           do k=ifrag_(1,jj),ifrag_(2,jj)
380             nhpb=nhpb+1
381             ihpb(nhpb)=j
382             jhpb(nhpb)=k
383             forcon(nhpb)=wpair_(i)
384             dhpb(nhpb)=dist(j,k)
385             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
386      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
387           enddo
388         enddo
389         endif
390       enddo 
391       do i=1,ndist_
392         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
393      &     ibecarb(i),forcon(nhpb+1)
394         if (forcon(nhpb+1).gt.0.0d0) then
395           nhpb=nhpb+1
396           if (ibecarb(i).gt.0) then
397             ihpb(i)=ihpb(i)+nres
398             jhpb(i)=jhpb(i)+nres
399           endif
400           if (dhpb(nhpb).eq.0.0d0) 
401      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
402         endif
403       enddo
404       do i=1,nhpb
405           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
406      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
407       enddo
408       call flush(iout)
409       return
410       end