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 if (nchain.gt.1) 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)
263 write (iout,*) "Calling init_int_table"
266 write (iout,'(/a,i3,a)') 'The chain contains',ns,
267 & ' disulfide-bridging cysteines.'
268 write (iout,'(20i4)') (iss(i),i=1,ns)
270 write(iout,*)"Running with dynamic disulfide-bond formation"
272 write (iout,'(/a/)') 'Pre-formed links are:'
278 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
279 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
280 & dhpb(i),ebr,forcon(i)
285 if (ns.gt.0.and.dyn_ss) then
289 forcon(i-nss)=forcon(i)
296 dyn_ss_mask(iss(i))=.true.
299 write (iout,*) "calling read_saxs_consrtr",nsaxs
300 if (nsaxs.gt.0) call read_saxs_constr
304 c-----------------------------------------------------------------------------
305 logical function seq_comp(itypea,itypeb,length)
307 integer length,itypea(length),itypeb(length)
310 if (itypea(i).ne.itypeb(i)) then
318 c-----------------------------------------------------------------------------
319 subroutine read_bridge
320 C Read information about disulfide bridges.
321 implicit real*8 (a-h,o-z)
323 include 'DIMENSIONS.ZSCOPT'
324 include 'COMMON.IOUNITS'
327 include 'COMMON.INTERACT'
328 include 'COMMON.NAMES'
329 include 'COMMON.CHAIN'
330 include 'COMMON.FFIELD'
331 include 'COMMON.SBRIDGE'
332 C Read bridging residues.
333 read (inp,*) ns,(iss(i),i=1,ns)
335 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
336 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
342 C Check whether the specified bridging residues are cystines.
344 if (itype(iss(i)).ne.1) then
345 write (iout,'(2a,i3,a)')
346 & 'Do you REALLY think that the residue ',
347 & restyp(itype(iss(i))),i,
348 & ' can form a disulfide bridge?!!!'
349 write (*,'(2a,i3,a)')
350 & 'Do you REALLY think that the residue ',
351 & restyp(itype(iss(i))),i,
352 & ' can form a disulfide bridge?!!!'
356 C Read preformed bridges.
358 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
359 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
362 C Check if the residues involved in bridges are in the specified list of
366 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
367 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
368 write (iout,'(a,i3,a)') 'Disulfide pair',i,
369 & ' contains residues present in other pairs.'
370 write (*,'(a,i3,a)') 'Disulfide pair',i,
371 & ' contains residues present in other pairs.'
376 if (ihpb(i).eq.iss(j)) goto 10
378 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
381 if (jhpb(i).eq.iss(j)) goto 20
383 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
396 c------------------------------------------------------------------------------
397 subroutine read_angles(kanal,iscor,energ,iprot,*)
398 implicit real*8 (a-h,o-z)
400 include 'DIMENSIONS.ZSCOPT'
401 include 'COMMON.INTERACT'
402 include 'COMMON.SBRIDGE'
405 include 'COMMON.CHAIN'
406 include 'COMMON.IOUNITS'
408 read(kanal,'(a80)',end=10,err=10) lineh
409 read(lineh(:5),*,err=8) ic
410 read(lineh(6:),*,err=8) energ
413 print *,'error, assuming e=1d10',lineh
417 read(lineh(18:),*,end=10,err=10) nss
419 read (lineh(20:),*,end=10,err=10)
420 & (IHPB(I),JHPB(I),I=1,NSS),iscor
422 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
423 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
426 c print *,"energy",energ," iscor",iscor
427 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
428 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
429 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
430 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
432 theta(i)=deg2rad*theta(i)
433 phi(i)=deg2rad*phi(i)
434 alph(i)=deg2rad*alph(i)
435 omeg(i)=deg2rad*omeg(i)
440 c-------------------------------------------------------------------------------
441 subroutine read_saxs_constr
442 implicit real*8 (a-h,o-z)
444 include 'DIMENSIONS.ZSCOPT'
445 include 'DIMENSIONS.FREE'
449 include 'COMMON.CONTROL'
450 include 'COMMON.CHAIN'
451 include 'COMMON.IOUNITS'
452 include 'COMMON.SBRIDGE'
453 include 'COMMON.SAXS'
454 double precision cm(3)
456 write (iout,*) "Calling read_saxs nsaxs",nsaxs
458 if (saxs_mode.eq.0) then
459 c SAXS distance distribution
461 read(inp,*) distsaxs(i),Psaxs(i)
465 Cnorm = Cnorm + Psaxs(i)
467 write (iout,*) "Cnorm",Cnorm
469 Psaxs(i)=Psaxs(i)/Cnorm
471 write (iout,*) "Normalized distance distribution from SAXS"
473 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
477 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
479 write (iout,*) "Wsaxs0",Wsaxs0
483 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
490 cm(j)=cm(j)+Csaxs(j,i)
498 Csaxs(j,i)=Csaxs(j,i)-cm(j)
501 write (iout,*) "SAXS sphere coordinates"
503 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)