5 implicit real*8 (a-h,o-z)
7 include 'DIMENSIONS.ZSCOPT'
8 include 'COMMON.IOUNITS'
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)
21 double precision x(maxvar)
22 character*320 controlcard,ucase
23 dimension itype_pdb(maxres)
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)
31 call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
32 call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
33 call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
34 call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
35 write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
36 & " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
37 r0_corr=cutoff_corr-delt_corr
38 call readi(controlcard,"NRES",nres,0)
39 iscode=index(controlcard,"ONE_LETTER")
41 write (iout,*) "Error: no residues in molecule"
44 if (nres.gt.maxres) then
45 write (iout,*) "Error: too many residues",nres,maxres
47 write(iout,*) 'nres=',nres
48 C Read sequence of the protein
50 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
52 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
54 C Convert sequence to numeric code
56 itype(i)=rescode(i,sequence(i),iscode)
58 write (iout,*) "Numeric code:"
59 write (iout,'(20i4)') (itype(i),i=1,nres)
62 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
64 if (itype(i).eq.21) then
68 else if (itype(i+1).ne.20) then
70 else if (itype(i).ne.20) then
79 if (with_dihed_constr) then
81 read (inp,*) ndih_constr
82 if (ndih_constr.gt.0) then
84 write (iout,*) 'FTORS',ftors
85 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
87 & 'There are',ndih_constr,' constraints on phi angles.'
89 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
92 phi0(i)=deg2rad*phi0(i)
93 drange(i)=deg2rad*drange(i)
101 if (itype(1).eq.21) nnt=2
102 if (itype(nres).eq.21) nct=nct-1
103 write(iout,*) 'NNT=',NNT,' NCT=',NCT
105 C Juyong:READ init_vars
106 C Initialize variables!
107 C Juyong:READ read_info
108 C READ fragment information!!
109 C both routines should be in dfa.F file!!
111 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
112 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
113 write (iout,*) "Calling init_dfa_vars"
116 write (iout,*) 'init_dfa_vars finished!'
119 write (iout,*) 'read_dfa_info finished!'
123 c Read distance restraints
124 if (constr_dist.gt.0) then
125 if (refstr) call read_ref_structure(*11)
126 call read_dist_constr
130 if (constr_homology.gt.0) then
131 call read_constr_homology
138 write (iout,'(/a,i3,a)') 'The chain contains',ns,
139 & ' disulfide-bridging cysteines.'
140 write (iout,'(20i4)') (iss(i),i=1,ns)
141 write (iout,'(/a/)') 'Pre-formed links are:'
147 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
148 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
149 & dhpb(i),ebr,forcon(i)
154 11 stop "Error reading reference structure"
156 c-----------------------------------------------------------------------------
157 logical function seq_comp(itypea,itypeb,length)
159 integer length,itypea(length),itypeb(length)
162 if (itypea(i).ne.itypeb(i)) then
170 c-----------------------------------------------------------------------------
171 subroutine read_bridge
172 C Read information about disulfide bridges.
173 implicit real*8 (a-h,o-z)
175 include 'DIMENSIONS.ZSCOPT'
176 include 'COMMON.IOUNITS'
179 include 'COMMON.INTERACT'
180 include 'COMMON.NAMES'
181 include 'COMMON.CHAIN'
182 include 'COMMON.FFIELD'
183 include 'COMMON.SBRIDGE'
184 C Read bridging residues.
185 read (inp,*) ns,(iss(i),i=1,ns)
187 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
188 C Check whether the specified bridging residues are cystines.
190 if (itype(iss(i)).ne.1) then
191 write (iout,'(2a,i3,a)')
192 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
193 & ' can form a disulfide bridge?!!!'
194 write (*,'(2a,i3,a)')
195 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
196 & ' can form a disulfide bridge?!!!'
200 C Read preformed bridges.
202 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
203 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
206 C Check if the residues involved in bridges are in the specified list of
210 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
211 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
212 write (iout,'(a,i3,a)') 'Disulfide pair',i,
213 & ' contains residues present in other pairs.'
214 write (*,'(a,i3,a)') 'Disulfide pair',i,
215 & ' contains residues present in other pairs.'
220 if (ihpb(i).eq.iss(j)) goto 10
222 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
225 if (jhpb(i).eq.iss(j)) goto 20
227 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
238 if (ns.gt.0.and.dyn_ss) then
239 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
242 C /06/28/2013 Adasko: nss number of full SS bonds
245 forcon(i-nss)=forcon(i)
252 dyn_ss_mask(iss(i))=.true.
253 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
254 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
259 c------------------------------------------------------------------------------
260 subroutine read_angles(kanal,iscor,energ,iprot,*)
261 implicit real*8 (a-h,o-z)
263 include 'DIMENSIONS.ZSCOPT'
264 include 'COMMON.INTERACT'
265 include 'COMMON.SBRIDGE'
268 include 'COMMON.CHAIN'
269 include 'COMMON.IOUNITS'
271 read(kanal,'(a80)',end=10,err=10) lineh
272 read(lineh(:5),*,err=8) ic
273 read(lineh(6:),*,err=8) energ
276 print *,'error, assuming e=1d10',lineh
280 read(lineh(18:),*,end=10,err=10) nss
282 read (lineh(20:),*,end=10,err=10)
283 & (IHPB(I),JHPB(I),I=1,NSS),iscor
285 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
286 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
289 c print *,"energy",energ," iscor",iscor
290 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
291 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
292 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
293 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
295 theta(i)=deg2rad*theta(i)
296 phi(i)=deg2rad*phi(i)
297 alph(i)=deg2rad*alph(i)
298 omeg(i)=deg2rad*omeg(i)
303 c-------------------------------------------------------------------------------
304 subroutine read_dist_constr
305 implicit real*8 (a-h,o-z)
307 include 'COMMON.CONTROL'
308 include 'COMMON.CHAIN'
309 include 'COMMON.IOUNITS'
310 include 'COMMON.SBRIDGE'
311 integer ifrag_(2,100),ipair_(2,100)
312 double precision wfrag_(100),wpair_(100)
313 character*500 controlcard
314 c write (iout,*) "Calling read_dist_constr"
315 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
317 call card_concat(controlcard,.true.)
318 call readi(controlcard,"NFRAG",nfrag_,0)
319 call readi(controlcard,"NPAIR",npair_,0)
320 call readi(controlcard,"NDIST",ndist_,0)
321 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
322 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
323 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
324 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
325 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
326 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
327 write (iout,*) "IFRAG"
329 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
331 write (iout,*) "IPAIR"
333 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
336 if (.not.refstr .and. nfrag_.gt.0) then
338 & "ERROR: no reference structure to compute distance restraints"
340 & "Restraints must be specified explicitly (NDIST=number)"
343 if (nfrag_.lt.2 .and. npair_.gt.0) then
344 write (iout,*) "ERROR: Less than 2 fragments specified",
345 & " but distance restraints between pairs requested"
350 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
351 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
352 & ifrag_(2,i)=nstart_sup+nsup-1
353 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
355 if (wfrag_(i).gt.0.0d0) then
356 do j=ifrag_(1,i),ifrag_(2,i)-1
358 write (iout,*) "j",j," k",k
360 if (constr_dist.eq.1) then
365 forcon(nhpb)=wfrag_(i)
366 else if (constr_dist.eq.2) then
367 if (ddjk.le.dist_cut) then
372 forcon(nhpb)=wfrag_(i)
379 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
381 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
382 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
388 if (wpair_(i).gt.0.0d0) then
396 do j=ifrag_(1,ii),ifrag_(2,ii)
397 do k=ifrag_(1,jj),ifrag_(2,jj)
401 forcon(nhpb)=wpair_(i)
403 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
404 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
410 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
411 & ibecarb(i),forcon(nhpb+1)
412 if (forcon(nhpb+1).gt.0.0d0) then
414 if (ibecarb(i).gt.0) then
418 if (dhpb(nhpb).eq.0.0d0)
419 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
423 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
424 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
432 c====-------------------------------------------------------------------
434 subroutine read_constr_homology
440 include 'COMMON.CONTROL'
441 include 'COMMON.CHAIN'
442 include 'COMMON.IOUNITS'
443 include 'COMMON.INTERACT'
445 double precision odl_temp,sigma_odl_temp
446 common /przechowalnia/ odl_temp(maxres,maxres,max_template),
447 & sigma_odl_temp(maxres,maxres,max_template)
449 character*24 model_ki_dist, model_ki_angle
450 character*500 controlcard
451 integer ki, i, j, k, l
452 logical lprn /.true./
454 call card_concat(controlcard,.true.)
455 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
456 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
457 write (iout,*) "wga_dist",waga_dist," waga_angle",waga_angle
463 do ki=1,constr_homology
464 sigma_odl_temp(i,j,ki)=0.0d0
465 odl_temp(i,j,ki)=0.0d0
470 do ki=1,constr_homology
472 sigma_dih(ki,i)=0.0d0
475 do ki=1,constr_homology
476 write(kic2,'(i2)') ki
477 if (ki.le.9) kic2="0"//kic2(2:2)
479 model_ki_dist="model"//kic2//".dist"
480 model_ki_angle="model"//kic2//".angle"
481 open (ientin,file=model_ki_dist,status='old')
482 do irec=1,maxdim !petla do czytania wiezow na odleglosc
483 read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
484 & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
485 odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
486 sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
487 & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
491 open (ientin,file=model_ki_angle,status='old')
492 do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
493 read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
494 & sigma_dih(ki,i+nnt-1)
495 if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
496 sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
502 write (iout,*) "nnt",nnt," nct",nct
506 c write (iout,*) "i",i," j",j," constr_homology",constr_homology
507 do while (ki.le.constr_homology .and.
508 & sigma_odl_temp(i,j,ki).le.0.0d0)
509 c write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
512 c write (iout,*) "ki",ki
513 if (ki.gt.constr_homology) cycle
517 do ki=1,constr_homology
518 odl(ki,ii)=odl_temp(i,j,ki)
519 sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
524 if (constr_homology.gt.0) call homology_partition
526 if (.not.lprn) return
527 write (iout,*) "Distance restraints from templates"
529 write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
530 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
532 write (iout,*) "Dihedral angle restraints from templates"
534 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
535 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
537 c write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
538 c write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)