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,
76 write(iout,*) "nres",nres," nchain",nchain
78 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
81 call chain_symmetry(nchain,nres,itype,chain_border,
82 & chain_length,npermchain,tabpermchain)
83 write(iout,*) "ireschain permutations"
85 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
88 write(iout,*) "residue permutations"
90 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
93 if (itype(1).eq.ntyp1) nnt=2
94 if (itype(nres).eq.ntyp1) nct=nct-1
95 write(iout,*) 'NNT=',NNT,' NCT=',NCT
96 if (with_dihed_constr) then
98 read (inp,*) ndih_constr
99 write (iout,*) "ndih_constr",ndih_constr
100 if (ndih_constr.gt.0) then
103 C write (iout,*) 'FTORS',ftors
104 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
107 & 'There are',ndih_constr,' restraints on gamma angles.'
109 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
113 phi0(i)=deg2rad*phi0(i)
114 drange(i)=deg2rad*drange(i)
116 else if (ndih_constr.lt.0) then
118 call card_concat(controlcard,.true.)
119 call reada(controlcard,"PHIHEL",phihel,50.0D0)
120 call reada(controlcard,"PHIBET",phibet,180.0D0)
121 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
122 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
123 call reada(controlcard,"WDIHC",wdihc,0.591d0)
124 write (iout,*) "Weight of the dihedral restraint term",wdihc
125 read(inp,'(9x,3f7.3)')
126 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
127 write (iout,*) "The secprob array"
129 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
133 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
134 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
135 ndih_constr=ndih_constr+1
136 idih_constr(ndih_constr)=i
139 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
140 sumv=sumv+vpsipred(j,ndih_constr)
143 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
145 phibound(1,ndih_constr)=phihel*deg2rad
146 phibound(2,ndih_constr)=phibet*deg2rad
147 sdihed(1,ndih_constr)=sigmahel*deg2rad
148 sdihed(2,ndih_constr)=sigmabet*deg2rad
152 & 'There are',ndih_constr,
153 & ' bimodal restraints on gamma angles.'
155 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
156 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
157 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
158 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
159 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
160 & (vpsipred(j,i),j=1,3)
166 if (with_theta_constr) then
167 C with_theta_constr is keyword allowing for occurance of theta constrains
168 read (inp,*) ntheta_constr
169 C ntheta_constr is the number of theta constrains
170 if (ntheta_constr.gt.0) then
172 read (inp,*) (itheta_constr(i),theta_constr0(i),
173 & theta_drange(i),for_thet_constr(i),
175 C the above code reads from 1 to ntheta_constr
176 C itheta_constr(i) residue i for which is theta_constr
177 C theta_constr0 the global minimum value
178 C theta_drange is range for which there is no energy penalty
179 C for_thet_constr is the force constant for quartic energy penalty
181 C if(me.eq.king.or..not.out1file)then
183 & 'There are',ntheta_constr,' constraints on phi angles.'
185 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
191 theta_constr0(i)=deg2rad*theta_constr0(i)
192 theta_drange(i)=deg2rad*theta_drange(i)
194 C if(me.eq.king.or..not.out1file)
195 C & write (iout,*) 'FTORS',ftors
196 C do i=1,ntheta_constr
197 C ii = itheta_constr(i)
198 C thetabound(1,ii) = phi0(i)-drange(i)
199 C thetabound(2,ii) = phi0(i)+drange(i)
201 endif ! ntheta_constr.gt.0
202 endif! with_theta_constr
203 if (constr_homology.gt.0) then
204 c write (iout,*) "About to call read_constr_homology"
206 call read_constr_homology
207 c write (iout,*) "Exit read_constr_homology"
209 if (indpdb.gt.0 .or. pdbref) then
213 cref(j,i)=crefjlee(j,i)
218 write (iout,*) "Array C"
220 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
221 & (c(j,i+nres),j=1,3)
223 write (iout,*) "Array Cref"
225 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
226 & (cref(j,i+nres),j=1,3)
230 call int_from_cart1(.false.)
231 call sc_loc_geom(.false.)
235 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
239 dc(j,i)=c(j,i+1)-c(j,i)
240 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
245 dc(j,i+nres)=c(j,i+nres)-c(j,i)
246 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
257 write (iout,'(/a,i3,a)') 'The chain contains',ns,
258 & ' disulfide-bridging cysteines.'
259 write (iout,'(20i4)') (iss(i),i=1,ns)
261 write(iout,*)"Running with dynamic disulfide-bond formation"
263 write (iout,'(/a/)') 'Pre-formed links are:'
269 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
270 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
271 & dhpb(i),ebr,forcon(i)
276 if (ns.gt.0.and.dyn_ss) then
280 forcon(i-nss)=forcon(i)
287 dyn_ss_mask(iss(i))=.true.
290 write (iout,*) "calling read_saxs_consrtr",nsaxs
291 if (nsaxs.gt.0) call read_saxs_constr
295 c-----------------------------------------------------------------------------
296 logical function seq_comp(itypea,itypeb,length)
298 integer length,itypea(length),itypeb(length)
301 if (itypea(i).ne.itypeb(i)) then
309 c-----------------------------------------------------------------------------
310 subroutine read_bridge
311 C Read information about disulfide bridges.
312 implicit real*8 (a-h,o-z)
314 include 'DIMENSIONS.ZSCOPT'
315 include 'COMMON.IOUNITS'
318 include 'COMMON.INTERACT'
319 include 'COMMON.NAMES'
320 include 'COMMON.CHAIN'
321 include 'COMMON.FFIELD'
322 include 'COMMON.SBRIDGE'
323 C Read bridging residues.
324 read (inp,*) ns,(iss(i),i=1,ns)
326 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
327 C Check whether the specified bridging residues are cystines.
329 if (itype(iss(i)).ne.1) then
330 write (iout,'(2a,i3,a)')
331 & 'Do you REALLY think that the residue ',
332 & restyp(itype(iss(i))),i,
333 & ' can form a disulfide bridge?!!!'
334 write (*,'(2a,i3,a)')
335 & 'Do you REALLY think that the residue ',
336 & restyp(itype(iss(i))),i,
337 & ' can form a disulfide bridge?!!!'
341 C Read preformed bridges.
343 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
344 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
347 C Check if the residues involved in bridges are in the specified list of
351 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
352 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
353 write (iout,'(a,i3,a)') 'Disulfide pair',i,
354 & ' contains residues present in other pairs.'
355 write (*,'(a,i3,a)') 'Disulfide pair',i,
356 & ' contains residues present in other pairs.'
361 if (ihpb(i).eq.iss(j)) goto 10
363 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
366 if (jhpb(i).eq.iss(j)) goto 20
368 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
381 c------------------------------------------------------------------------------
382 subroutine read_angles(kanal,iscor,energ,iprot,*)
383 implicit real*8 (a-h,o-z)
385 include 'DIMENSIONS.ZSCOPT'
386 include 'COMMON.INTERACT'
387 include 'COMMON.SBRIDGE'
390 include 'COMMON.CHAIN'
391 include 'COMMON.IOUNITS'
393 read(kanal,'(a80)',end=10,err=10) lineh
394 read(lineh(:5),*,err=8) ic
395 read(lineh(6:),*,err=8) energ
398 print *,'error, assuming e=1d10',lineh
402 read(lineh(18:),*,end=10,err=10) nss
404 read (lineh(20:),*,end=10,err=10)
405 & (IHPB(I),JHPB(I),I=1,NSS),iscor
407 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
408 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
411 c print *,"energy",energ," iscor",iscor
412 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
413 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
414 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
415 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
417 theta(i)=deg2rad*theta(i)
418 phi(i)=deg2rad*phi(i)
419 alph(i)=deg2rad*alph(i)
420 omeg(i)=deg2rad*omeg(i)
425 c-------------------------------------------------------------------------------
426 subroutine read_saxs_constr
427 implicit real*8 (a-h,o-z)
429 include 'DIMENSIONS.ZSCOPT'
430 include 'DIMENSIONS.FREE'
434 include 'COMMON.CONTROL'
435 include 'COMMON.CHAIN'
436 include 'COMMON.IOUNITS'
437 include 'COMMON.SBRIDGE'
438 include 'COMMON.SAXS'
439 double precision cm(3)
441 write (iout,*) "Calling read_saxs nsaxs",nsaxs
443 if (saxs_mode.eq.0) then
444 c SAXS distance distribution
446 read(inp,*) distsaxs(i),Psaxs(i)
450 Cnorm = Cnorm + Psaxs(i)
452 write (iout,*) "Cnorm",Cnorm
454 Psaxs(i)=Psaxs(i)/Cnorm
456 write (iout,*) "Normalized distance distribution from SAXS"
458 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
462 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
464 write (iout,*) "Wsaxs0",Wsaxs0
468 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
475 cm(j)=cm(j)+Csaxs(j,i)
483 Csaxs(j,i)=Csaxs(j,i)-cm(j)
486 write (iout,*) "SAXS sphere coordinates"
488 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)