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 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")
35 write (iout,*) "Error: no residues in molecule"
38 if (nres.gt.maxres) then
39 write (iout,*) "Error: too many residues",nres,maxres
41 write(iout,*) 'nres=',nres
42 C Read sequence of the protein
44 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
46 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
48 C Convert sequence to numeric code
50 itype(i)=rescode(i,sequence(i),iscode)
52 write (iout,*) "Numeric code:"
53 write (iout,'(20i4)') (itype(i),i=1,nres)
56 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
58 if (itype(i).eq.ntyp1) then
62 else if (iabs(itype(i+1)).ne.20) then
64 else if (iabs(itype(i)).ne.20) then
73 write (iout,*) i,itype(i),itel(i)
78 if (itype(1).eq.ntyp1) nnt=2
79 if (itype(nres).eq.ntyp1) nct=nct-1
80 write(iout,*) 'NNT=',NNT,' NCT=',NCT
82 if (with_dihed_constr) then
84 read (inp,*) ndih_constr
85 write (iout,*) "ndih_constr",ndih_constr
86 if (ndih_constr.gt.0) then
89 C write (iout,*) 'FTORS',ftors
90 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
93 & 'There are',ndih_constr,' restraints on gamma angles.'
95 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
99 phi0(i)=deg2rad*phi0(i)
100 drange(i)=deg2rad*drange(i)
102 else if (ndih_constr.lt.0) then
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"
115 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
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
125 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
126 sumv=sumv+vpsipred(j,ndih_constr)
129 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
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
138 & 'There are',ndih_constr,
139 & ' bimodal restraints on gamma angles.'
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)
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
158 read (inp,*) (itheta_constr(i),theta_constr0(i),
159 & theta_drange(i),for_thet_constr(i),
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
167 C if(me.eq.king.or..not.out1file)then
169 & 'There are',ntheta_constr,' constraints on phi angles.'
171 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
177 theta_constr0(i)=deg2rad*theta_constr0(i)
178 theta_drange(i)=deg2rad*theta_drange(i)
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)
187 endif ! ntheta_constr.gt.0
188 endif! with_theta_constr
192 write (iout,'(/a,i3,a)') 'The chain contains',ns,
193 & ' disulfide-bridging cysteines.'
194 write (iout,'(20i4)') (iss(i),i=1,ns)
196 write(iout,*)"Running with dynamic disulfide-bond formation"
198 write (iout,'(/a/)') 'Pre-formed links are:'
204 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
205 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
206 & dhpb(i),ebr,forcon(i)
211 if (ns.gt.0.and.dyn_ss) then
215 forcon(i-nss)=forcon(i)
222 dyn_ss_mask(iss(i))=.true.
225 write (iout,*) "calling read_saxs_consrtr",nsaxs
226 if (nsaxs.gt.0) call read_saxs_constr
230 c-----------------------------------------------------------------------------
231 logical function seq_comp(itypea,itypeb,length)
233 integer length,itypea(length),itypeb(length)
236 if (itypea(i).ne.itypeb(i)) then
244 c-----------------------------------------------------------------------------
245 subroutine read_bridge
246 C Read information about disulfide bridges.
247 implicit real*8 (a-h,o-z)
249 include 'DIMENSIONS.ZSCOPT'
250 include 'COMMON.IOUNITS'
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)
261 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
262 C Check whether the specified bridging residues are cystines.
264 if (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?!!!'
276 C Read preformed bridges.
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)
282 C Check if the residues involved in bridges are in the specified list of
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.'
296 if (ihpb(i).eq.iss(j)) goto 10
298 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
301 if (jhpb(i).eq.iss(j)) goto 20
303 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
316 c------------------------------------------------------------------------------
317 subroutine read_angles(kanal,iscor,energ,iprot,*)
318 implicit real*8 (a-h,o-z)
320 include 'DIMENSIONS.ZSCOPT'
321 include 'COMMON.INTERACT'
322 include 'COMMON.SBRIDGE'
325 include 'COMMON.CHAIN'
326 include 'COMMON.IOUNITS'
328 read(kanal,'(a80)',end=10,err=10) lineh
329 read(lineh(:5),*,err=8) ic
330 read(lineh(6:),*,err=8) energ
333 print *,'error, assuming e=1d10',lineh
337 read(lineh(18:),*,end=10,err=10) nss
339 read (lineh(20:),*,end=10,err=10)
340 & (IHPB(I),JHPB(I),I=1,NSS),iscor
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),
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)
352 theta(i)=deg2rad*theta(i)
353 phi(i)=deg2rad*phi(i)
354 alph(i)=deg2rad*alph(i)
355 omeg(i)=deg2rad*omeg(i)
360 c-------------------------------------------------------------------------------
361 subroutine read_saxs_constr
362 implicit real*8 (a-h,o-z)
364 include 'DIMENSIONS.ZSCOPT'
365 include 'DIMENSIONS.FREE'
369 include 'COMMON.CONTROL'
370 include 'COMMON.CHAIN'
371 include 'COMMON.IOUNITS'
372 include 'COMMON.SBRIDGE'
373 double precision cm(3)
375 write (iout,*) "Calling read_saxs nsaxs",nsaxs
377 if (saxs_mode.eq.0) then
378 c SAXS distance distribution
380 read(inp,*) distsaxs(i),Psaxs(i)
384 Cnorm = Cnorm + Psaxs(i)
386 write (iout,*) "Cnorm",Cnorm
388 Psaxs(i)=Psaxs(i)/Cnorm
390 write (iout,*) "Normalized distance distribution from SAXS"
392 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
396 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
398 write (iout,*) "Wsaxs0",Wsaxs0
402 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
409 cm(j)=cm(j)+Csaxs(j,i)
417 Csaxs(j,i)=Csaxs(j,i)-cm(j)
420 write (iout,*) "SAXS sphere coordinates"
422 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)