update
[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       double precision secprob(3,maxdih_constr),phihel,phibet
26       call card_concat(controlcard,.true.)
27       call reada(controlcard,'SCAL14',scal14,0.4d0)
28       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
29       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
30       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
31       r0_corr=cutoff_corr-delt_corr
32       call readi(controlcard,"NRES",nres,0)
33       iscode=index(controlcard,"ONE_LETTER")
34       if (nres.le.0) then
35         write (iout,*) "Error: no residues in molecule"
36         return1
37       endif
38       if (nres.gt.maxres) then
39         write (iout,*) "Error: too many residues",nres,maxres
40       endif
41       write(iout,*) 'nres=',nres
42 C Read sequence of the protein
43       if (iscode.gt.0) then
44         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
45       else
46         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
47       endif
48 C Convert sequence to numeric code
49       do i=1,nres
50         itype(i)=rescode(i,sequence(i),iscode)
51       enddo
52       write (iout,*) "Numeric code:"
53       write (iout,'(20i4)') (itype(i),i=1,nres)
54       do i=1,nres-1
55 #ifdef PROCOR
56         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
57 #else
58         if (itype(i).eq.ntyp1) then
59 #endif
60           itel(i)=0
61 #ifdef PROCOR
62         else if (iabs(itype(i+1)).ne.20) then
63 #else
64         else if (iabs(itype(i)).ne.20) then
65 #endif
66           itel(i)=1
67         else
68           itel(i)=2
69         endif
70       enddo
71        write (iout,*) "ITEL"
72        do i=1,nres-1
73          write (iout,*) i,itype(i),itel(i)
74        enddo
75       call read_bridge
76       nnt=1
77       nct=nres
78       if (itype(1).eq.ntyp1) nnt=2
79       if (itype(nres).eq.ntyp1) nct=nct-1
80       write(iout,*) 'NNT=',NNT,' NCT=',NCT
81
82       if (with_dihed_constr) then
83
84       read (inp,*) ndih_constr
85       write (iout,*) "ndih_constr",ndih_constr
86       if (ndih_constr.gt.0) then
87          raw_psipred=.false.
88 C        read (inp,*) ftors
89 C        write (iout,*) 'FTORS',ftors
90         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
91      &   i=1,ndih_constr)
92         write (iout,*)
93      &   'There are',ndih_constr,' restraints on gamma angles.'
94         do i=1,ndih_constr
95           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
96      &    ftors(i)
97         enddo
98         do i=1,ndih_constr
99           phi0(i)=deg2rad*phi0(i)
100           drange(i)=deg2rad*drange(i)
101         enddo
102       else if (ndih_constr.lt.0) then
103         raw_psipred=.true.
104         call card_concat(controlcard,.true.)
105         call reada(controlcard,"PHIHEL",phihel,50.0D0)
106         call reada(controlcard,"PHIBET",phibet,180.0D0)
107         call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
108         call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
109         call reada(controlcard,"WDIHC",wdihc,0.591d0)
110         write (iout,*) "Weight of the dihedral restraint term",wdihc
111         read(inp,'(9x,3f7.3)')
112      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
113         write (iout,*) "The secprob array"
114         do i=nnt,nct
115           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
116         enddo
117         ndih_constr=0
118         do i=nnt+3,nct
119           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
120      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
121             ndih_constr=ndih_constr+1
122             idih_constr(ndih_constr)=i
123             sumv=0.0d0
124             do j=1,3
125               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
126               sumv=sumv+vpsipred(j,ndih_constr)
127             enddo
128             do j=1,3
129               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
130             enddo
131             phibound(1,ndih_constr)=phihel*deg2rad
132             phibound(2,ndih_constr)=phibet*deg2rad
133             sdihed(1,ndih_constr)=sigmahel*deg2rad
134             sdihed(2,ndih_constr)=sigmabet*deg2rad
135           endif
136         enddo
137         write (iout,*)
138      &   'There are',ndih_constr,
139      &   ' bimodal restraints on gamma angles.'
140         do i=1,ndih_constr
141           write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
142      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
143      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
144      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
145      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
146      &      (vpsipred(j,i),j=1,3)
147         enddo
148
149       endif
150
151       endif
152       if (with_theta_constr) then
153 C with_theta_constr is keyword allowing for occurance of theta constrains
154       read (inp,*) ntheta_constr
155 C ntheta_constr is the number of theta constrains
156       if (ntheta_constr.gt.0) then
157 C        read (inp,*) ftors
158         read (inp,*) (itheta_constr(i),theta_constr0(i),
159      &  theta_drange(i),for_thet_constr(i),
160      &  i=1,ntheta_constr)
161 C the above code reads from 1 to ntheta_constr 
162 C itheta_constr(i) residue i for which is theta_constr
163 C theta_constr0 the global minimum value
164 C theta_drange is range for which there is no energy penalty
165 C for_thet_constr is the force constant for quartic energy penalty
166 C E=k*x**4 
167 C        if(me.eq.king.or..not.out1file)then
168          write (iout,*)
169      &   'There are',ntheta_constr,' constraints on phi angles.'
170          do i=1,ntheta_constr
171           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
172      &    theta_drange(i),
173      &    for_thet_constr(i)
174          enddo
175 C        endif
176         do i=1,ntheta_constr
177           theta_constr0(i)=deg2rad*theta_constr0(i)
178           theta_drange(i)=deg2rad*theta_drange(i)
179         enddo
180 C        if(me.eq.king.or..not.out1file)
181 C     &   write (iout,*) 'FTORS',ftors
182 C        do i=1,ntheta_constr
183 C          ii = itheta_constr(i)
184 C          thetabound(1,ii) = phi0(i)-drange(i)
185 C          thetabound(2,ii) = phi0(i)+drange(i)
186 C        enddo
187       endif ! ntheta_constr.gt.0
188       endif! with_theta_constr
189       call setup_var
190       call init_int_table
191       if (ns.gt.0) then
192         write (iout,'(/a,i3,a)') 'The chain contains',ns,
193      &  ' disulfide-bridging cysteines.'
194         write (iout,'(20i4)') (iss(i),i=1,ns)
195        if (dyn_ss) then
196           write(iout,*)"Running with dynamic disulfide-bond formation"
197        else
198         write (iout,'(/a/)') 'Pre-formed links are:' 
199         do i=1,nss
200           i1=ihpb(i)-nres
201           i2=jhpb(i)-nres
202           it1=itype(i1)
203           it2=itype(i2)
204          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
205      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
206      &    dhpb(i),ebr,forcon(i)
207         enddo
208       endif
209       endif
210       write (iout,'(a)')
211       write (iout,*) "setting ss mask dyn_ss",dyn_ss
212       if (ns.gt.0.and.dyn_ss) then
213           do i=nss+1,nhpb
214             ihpb(i-nss)=ihpb(i)
215             jhpb(i-nss)=jhpb(i)
216             forcon(i-nss)=forcon(i)
217             dhpb(i-nss)=dhpb(i)
218           enddo
219           nhpb=nhpb-nss
220           nss=0
221           call hpb_partition
222           do i=1,ns
223             dyn_ss_mask(iss(i))=.true.
224             write (iout,*) "i",i," iss",iss(i),
225      &        " mask",dyn_ss_mask(iss(i))
226           enddo
227       endif
228       return
229       end
230 c-----------------------------------------------------------------------------
231       logical function seq_comp(itypea,itypeb,length)
232       implicit none
233       integer length,itypea(length),itypeb(length)
234       integer i
235       do i=1,length
236         if (itypea(i).ne.itypeb(i)) then
237           seq_comp=.false.
238           return
239         endif
240       enddo
241       seq_comp=.true.
242       return
243       end
244 c-----------------------------------------------------------------------------
245       subroutine read_bridge
246 C Read information about disulfide bridges.
247       implicit real*8 (a-h,o-z)
248       include 'DIMENSIONS'
249       include 'DIMENSIONS.ZSCOPT'
250       include 'COMMON.IOUNITS'
251       include 'COMMON.GEO'
252       include 'COMMON.VAR'
253       include 'COMMON.INTERACT'
254       include 'COMMON.NAMES'
255       include 'COMMON.CHAIN'
256       include 'COMMON.FFIELD'
257       include 'COMMON.SBRIDGE'
258 C Read bridging residues.
259       read (inp,*) ns,(iss(i),i=1,ns)
260       print *,'ns=',ns
261       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
262 C Check whether the specified bridging residues are cystines.
263       do i=1,ns
264         if (iabs(itype(iss(i))).ne.1) then
265           write (iout,'(2a,i3,a)') 
266      &   'Do you REALLY think that the residue ',
267      &    restyp(itype(iss(i))),i,
268      &   ' can form a disulfide bridge?!!!'
269           write (*,'(2a,i3,a)') 
270      &   'Do you REALLY think that the residue ',
271      &    restyp(itype(iss(i))),i,
272      &   ' can form a disulfide bridge?!!!'
273          stop
274         endif
275       enddo
276 C Read preformed bridges.
277       if (ns.gt.0) then
278       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
279       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
280       if (nss.gt.0) then
281         nhpb=nss
282 C Check if the residues involved in bridges are in the specified list of
283 C bridging residues.
284         do i=1,nss
285           do j=1,i-1
286             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
287      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
288               write (iout,'(a,i3,a)') 'Disulfide pair',i,
289      &      ' contains residues present in other pairs.'
290               write (*,'(a,i3,a)') 'Disulfide pair',i,
291      &      ' contains residues present in other pairs.'
292               stop 
293             endif
294           enddo
295           do j=1,ns
296             if (ihpb(i).eq.iss(j)) goto 10
297           enddo
298           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
299    10     continue
300           do j=1,ns
301             if (jhpb(i).eq.iss(j)) goto 20
302           enddo
303           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
304    20     continue
305           dhpb(i)=dbr
306           forcon(i)=fbr
307         enddo
308         do i=1,nss
309           ihpb(i)=ihpb(i)+nres
310           jhpb(i)=jhpb(i)+nres
311         enddo
312       endif
313       endif
314       return
315       end
316 c------------------------------------------------------------------------------
317       subroutine read_angles(kanal,iscor,energ,iprot,*)
318       implicit real*8 (a-h,o-z)
319       include 'DIMENSIONS'
320       include 'DIMENSIONS.ZSCOPT'
321       include 'COMMON.INTERACT'
322       include 'COMMON.SBRIDGE'
323       include 'COMMON.GEO'
324       include 'COMMON.VAR'
325       include 'COMMON.CHAIN'
326       include 'COMMON.IOUNITS'
327       character*80 lineh
328       read(kanal,'(a80)',end=10,err=10) lineh
329       read(lineh(:5),*,err=8) ic
330       read(lineh(6:),*,err=8) energ
331       goto 9
332     8 ic=1
333       print *,'error, assuming e=1d10',lineh
334       energ=1d10
335       nss=0
336     9 continue
337       read(lineh(18:),*,end=10,err=10) nss
338       IF (NSS.LT.9) THEN
339         read (lineh(20:),*,end=10,err=10)
340      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
341       ELSE
342         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
343         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
344      &    I=9,NSS),iscor
345       ENDIF
346 c      print *,"energy",energ," iscor",iscor
347       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
348       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
349       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
350       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
351       do i=1,nres
352         theta(i)=deg2rad*theta(i)
353         phi(i)=deg2rad*phi(i)
354         alph(i)=deg2rad*alph(i)
355         omeg(i)=deg2rad*omeg(i)
356       enddo
357       return
358    10 return1
359       end