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'
20 character*4 sequence(maxres)
22 double precision x(maxvar)
23 character*320 controlcard,ucase
24 dimension itype_pdb(maxres)
26 double precision secprob(3,maxdih_constr),phihel,phibet
27 call card_concat(controlcard,.true.)
28 call readi(controlcard,"NRES",nres,0)
29 iscode=index(controlcard,"ONE_LETTER")
31 write (iout,*) "Error: no residues in molecule"
34 if (nres.gt.maxres) then
35 write (iout,*) "Error: too many residues",nres,maxres
37 write(iout,*) 'nres=',nres
38 C Read sequence of the protein
40 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
42 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
44 C Convert sequence to numeric code
46 itype(i)=rescode(i,sequence(i),iscode)
48 write (iout,*) "Numeric code:"
49 write (iout,'(20i4)') (itype(i),i=1,nres)
52 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
54 if (itype(i).eq.ntyp1) then
58 else if (iabs(itype(i+1)).ne.20) then
60 else if (iabs(itype(i)).ne.20) then
69 write (iout,*) i,itype(i),itel(i)
74 call seq2chains(nres,itype,nchain,chain_length,chain_border,
77 chain_border1(2,1)=chain_border(2,1)+1
79 chain_border1(1,i)=chain_border(1,i)-1
80 chain_border1(2,i)=chain_border(2,i)+1
82 chain_border1(1,nchain)=chain_border(1,nchain)-1
83 chain_border1(2,nchain)=nres
84 write(iout,*) "nres",nres," nchain",nchain
86 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
87 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
89 call chain_symmetry(nchain,nres,itype,chain_border,
90 & chain_length,npermchain,tabpermchain)
91 write(iout,*) "ireschain permutations"
93 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
96 write(iout,*) "residue permutations"
98 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
101 if (itype(1).eq.ntyp1) nnt=2
102 if (itype(nres).eq.ntyp1) nct=nct-1
103 write(iout,*) 'NNT=',NNT,' NCT=',NCT
104 if (with_dihed_constr) then
106 read (inp,*) ndih_constr
107 write (iout,*) "ndih_constr",ndih_constr
108 if (ndih_constr.gt.0) then
111 C write (iout,*) 'FTORS',ftors
112 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
115 & 'There are',ndih_constr,' restraints on gamma angles.'
117 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
121 phi0(i)=deg2rad*phi0(i)
122 drange(i)=deg2rad*drange(i)
124 else if (ndih_constr.lt.0) then
126 call card_concat(controlcard,.true.)
127 call reada(controlcard,"PHIHEL",phihel,50.0D0)
128 call reada(controlcard,"PHIBET",phibet,180.0D0)
129 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
130 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
131 call reada(controlcard,"WDIHC",wdihc,0.591d0)
132 write (iout,*) "Weight of the dihedral restraint term",wdihc
133 read(inp,'(9x,3f7.3)')
134 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
135 write (iout,*) "The secprob array"
137 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
141 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
142 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
143 ndih_constr=ndih_constr+1
144 idih_constr(ndih_constr)=i
147 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
148 sumv=sumv+vpsipred(j,ndih_constr)
151 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
153 phibound(1,ndih_constr)=phihel*deg2rad
154 phibound(2,ndih_constr)=phibet*deg2rad
155 sdihed(1,ndih_constr)=sigmahel*deg2rad
156 sdihed(2,ndih_constr)=sigmabet*deg2rad
160 & 'There are',ndih_constr,
161 & ' bimodal restraints on gamma angles.'
163 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
164 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
165 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
166 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
167 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
168 & (vpsipred(j,i),j=1,3)
174 if (with_theta_constr) then
175 C with_theta_constr is keyword allowing for occurance of theta constrains
176 read (inp,*) ntheta_constr
177 C ntheta_constr is the number of theta constrains
178 if (ntheta_constr.gt.0) then
180 read (inp,*) (itheta_constr(i),theta_constr0(i),
181 & theta_drange(i),for_thet_constr(i),
183 C the above code reads from 1 to ntheta_constr
184 C itheta_constr(i) residue i for which is theta_constr
185 C theta_constr0 the global minimum value
186 C theta_drange is range for which there is no energy penalty
187 C for_thet_constr is the force constant for quartic energy penalty
189 C if(me.eq.king.or..not.out1file)then
191 & 'There are',ntheta_constr,' constraints on phi angles.'
193 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
199 theta_constr0(i)=deg2rad*theta_constr0(i)
200 theta_drange(i)=deg2rad*theta_drange(i)
202 C if(me.eq.king.or..not.out1file)
203 C & write (iout,*) 'FTORS',ftors
204 C do i=1,ntheta_constr
205 C ii = itheta_constr(i)
206 C thetabound(1,ii) = phi0(i)-drange(i)
207 C thetabound(2,ii) = phi0(i)+drange(i)
209 endif ! ntheta_constr.gt.0
210 endif! with_theta_constr
211 if (constr_homology.gt.0) then
212 c write (iout,*) "About to call read_constr_homology"
214 call read_constr_homology
215 c write (iout,*) "Exit read_constr_homology"
217 if (indpdb.gt.0 .or. pdbref) then
221 cref(j,i)=crefjlee(j,i)
226 write (iout,*) "Array C"
228 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
229 & (c(j,i+nres),j=1,3)
231 write (iout,*) "Array Cref"
233 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
234 & (cref(j,i+nres),j=1,3)
238 call int_from_cart1(.false.)
239 call sc_loc_geom(.false.)
243 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
247 dc(j,i)=c(j,i+1)-c(j,i)
248 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
253 dc(j,i+nres)=c(j,i+nres)-c(j,i)
254 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
265 write (iout,'(/a,i3,a)') 'The chain contains',ns,
266 & ' disulfide-bridging cysteines.'
267 write (iout,'(20i4)') (iss(i),i=1,ns)
269 write(iout,*)"Running with dynamic disulfide-bond formation"
271 write (iout,'(/a/)') 'Pre-formed links are:'
277 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
278 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
279 & dhpb(i),ebr,forcon(i)
284 if (ns.gt.0.and.dyn_ss) then
288 forcon(i-nss)=forcon(i)
295 dyn_ss_mask(iss(i))=.true.
298 write (iout,*) "calling read_saxs_consrtr",nsaxs
299 if (nsaxs.gt.0) call read_saxs_constr
303 c-----------------------------------------------------------------------------
304 logical function seq_comp(itypea,itypeb,length)
306 integer length,itypea(length),itypeb(length)
309 if (itypea(i).ne.itypeb(i)) then
317 c-----------------------------------------------------------------------------
318 subroutine read_bridge
319 C Read information about disulfide bridges.
320 implicit real*8 (a-h,o-z)
322 include 'DIMENSIONS.ZSCOPT'
323 include 'COMMON.IOUNITS'
326 include 'COMMON.INTERACT'
327 include 'COMMON.NAMES'
328 include 'COMMON.CHAIN'
329 include 'COMMON.FFIELD'
330 include 'COMMON.SBRIDGE'
331 C Read bridging residues.
332 read (inp,*) ns,(iss(i),i=1,ns)
334 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
335 C Check whether the specified bridging residues are cystines.
337 if (itype(iss(i)).ne.1) then
338 write (iout,'(2a,i3,a)')
339 & 'Do you REALLY think that the residue ',
340 & restyp(itype(iss(i))),i,
341 & ' can form a disulfide bridge?!!!'
342 write (*,'(2a,i3,a)')
343 & 'Do you REALLY think that the residue ',
344 & restyp(itype(iss(i))),i,
345 & ' can form a disulfide bridge?!!!'
349 C Read preformed bridges.
351 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
352 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
355 C Check if the residues involved in bridges are in the specified list of
359 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
360 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
361 write (iout,'(a,i3,a)') 'Disulfide pair',i,
362 & ' contains residues present in other pairs.'
363 write (*,'(a,i3,a)') 'Disulfide pair',i,
364 & ' contains residues present in other pairs.'
369 if (ihpb(i).eq.iss(j)) goto 10
371 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
374 if (jhpb(i).eq.iss(j)) goto 20
376 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
389 c------------------------------------------------------------------------------
390 subroutine read_angles(kanal,iscor,energ,iprot,*)
391 implicit real*8 (a-h,o-z)
393 include 'DIMENSIONS.ZSCOPT'
394 include 'COMMON.INTERACT'
395 include 'COMMON.SBRIDGE'
398 include 'COMMON.CHAIN'
399 include 'COMMON.IOUNITS'
401 read(kanal,'(a80)',end=10,err=10) lineh
402 read(lineh(:5),*,err=8) ic
403 read(lineh(6:),*,err=8) energ
406 print *,'error, assuming e=1d10',lineh
410 read(lineh(18:),*,end=10,err=10) nss
412 read (lineh(20:),*,end=10,err=10)
413 & (IHPB(I),JHPB(I),I=1,NSS),iscor
415 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
416 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
419 c print *,"energy",energ," iscor",iscor
420 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
421 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
422 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
423 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
425 theta(i)=deg2rad*theta(i)
426 phi(i)=deg2rad*phi(i)
427 alph(i)=deg2rad*alph(i)
428 omeg(i)=deg2rad*omeg(i)
433 c-------------------------------------------------------------------------------
434 subroutine read_saxs_constr
435 implicit real*8 (a-h,o-z)
437 include 'DIMENSIONS.ZSCOPT'
438 include 'DIMENSIONS.FREE'
442 include 'COMMON.CONTROL'
443 include 'COMMON.CHAIN'
444 include 'COMMON.IOUNITS'
445 include 'COMMON.SBRIDGE'
446 include 'COMMON.SAXS'
447 double precision cm(3)
449 write (iout,*) "Calling read_saxs nsaxs",nsaxs
451 if (saxs_mode.eq.0) then
452 c SAXS distance distribution
454 read(inp,*) distsaxs(i),Psaxs(i)
458 Cnorm = Cnorm + Psaxs(i)
460 write (iout,*) "Cnorm",Cnorm
462 Psaxs(i)=Psaxs(i)/Cnorm
464 write (iout,*) "Normalized distance distribution from SAXS"
466 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
470 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
472 write (iout,*) "Wsaxs0",Wsaxs0
476 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
483 cm(j)=cm(j)+Csaxs(j,i)
491 Csaxs(j,i)=Csaxs(j,i)-cm(j)
494 write (iout,*) "SAXS sphere coordinates"
496 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)