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 call seq2chains(nres,itype,nchain,chain_length,chain_border,
80 write(iout,*) "nres",nres," nchain",nchain
82 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
85 call chain_symmetry(nchain,nres,itype,chain_border,
86 & chain_length,npermchain,tabpermchain)
87 write(iout,*) "ireschain permutations"
89 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
92 write(iout,*) "residue permutations"
94 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
97 if (itype(1).eq.ntyp1) nnt=2
98 if (itype(nres).eq.ntyp1) nct=nct-1
99 write(iout,*) 'NNT=',NNT,' NCT=',NCT
101 if (with_dihed_constr) then
103 read (inp,*) ndih_constr
104 write (iout,*) "ndih_constr",ndih_constr
105 if (ndih_constr.gt.0) then
108 C write (iout,*) 'FTORS',ftors
109 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
112 & 'There are',ndih_constr,' restraints on gamma angles.'
114 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
118 phi0(i)=deg2rad*phi0(i)
119 drange(i)=deg2rad*drange(i)
121 else if (ndih_constr.lt.0) then
123 call card_concat(controlcard,.true.)
124 call reada(controlcard,"PHIHEL",phihel,50.0D0)
125 call reada(controlcard,"PHIBET",phibet,180.0D0)
126 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
127 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
128 call reada(controlcard,"WDIHC",wdihc,0.591d0)
129 write (iout,*) "Weight of the dihedral restraint term",wdihc
130 read(inp,'(9x,3f7.3)')
131 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
132 write (iout,*) "The secprob array"
134 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
138 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
139 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
140 ndih_constr=ndih_constr+1
141 idih_constr(ndih_constr)=i
144 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
145 sumv=sumv+vpsipred(j,ndih_constr)
148 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
150 phibound(1,ndih_constr)=phihel*deg2rad
151 phibound(2,ndih_constr)=phibet*deg2rad
152 sdihed(1,ndih_constr)=sigmahel*deg2rad
153 sdihed(2,ndih_constr)=sigmabet*deg2rad
157 & 'There are',ndih_constr,
158 & ' bimodal restraints on gamma angles.'
160 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
161 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
162 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
163 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
164 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
165 & (vpsipred(j,i),j=1,3)
171 if (with_theta_constr) then
172 C with_theta_constr is keyword allowing for occurance of theta constrains
173 read (inp,*) ntheta_constr
174 C ntheta_constr is the number of theta constrains
175 if (ntheta_constr.gt.0) then
177 read (inp,*) (itheta_constr(i),theta_constr0(i),
178 & theta_drange(i),for_thet_constr(i),
180 C the above code reads from 1 to ntheta_constr
181 C itheta_constr(i) residue i for which is theta_constr
182 C theta_constr0 the global minimum value
183 C theta_drange is range for which there is no energy penalty
184 C for_thet_constr is the force constant for quartic energy penalty
186 C if(me.eq.king.or..not.out1file)then
188 & 'There are',ntheta_constr,' constraints on phi angles.'
190 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
196 theta_constr0(i)=deg2rad*theta_constr0(i)
197 theta_drange(i)=deg2rad*theta_drange(i)
199 C if(me.eq.king.or..not.out1file)
200 C & write (iout,*) 'FTORS',ftors
201 C do i=1,ntheta_constr
202 C ii = itheta_constr(i)
203 C thetabound(1,ii) = phi0(i)-drange(i)
204 C thetabound(2,ii) = phi0(i)+drange(i)
206 endif ! ntheta_constr.gt.0
207 endif! with_theta_constr
211 write (iout,'(/a,i3,a)') 'The chain contains',ns,
212 & ' disulfide-bridging cysteines.'
213 write (iout,'(20i4)') (iss(i),i=1,ns)
215 write(iout,*)"Running with dynamic disulfide-bond formation"
217 write (iout,'(/a/)') 'Pre-formed links are:'
223 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
224 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
225 & dhpb(i),ebr,forcon(i)
230 if (ns.gt.0.and.dyn_ss) then
234 forcon(i-nss)=forcon(i)
241 dyn_ss_mask(iss(i))=.true.
244 write (iout,*) "calling read_saxs_consrtr",nsaxs
245 if (nsaxs.gt.0) call read_saxs_constr
249 c-----------------------------------------------------------------------------
250 logical function seq_comp(itypea,itypeb,length)
252 integer length,itypea(length),itypeb(length)
255 if (itypea(i).ne.itypeb(i)) then
263 c-----------------------------------------------------------------------------
264 subroutine read_bridge
265 C Read information about disulfide bridges.
266 implicit real*8 (a-h,o-z)
268 include 'DIMENSIONS.ZSCOPT'
269 include 'COMMON.IOUNITS'
272 include 'COMMON.INTERACT'
273 include 'COMMON.NAMES'
274 include 'COMMON.CHAIN'
275 include 'COMMON.FFIELD'
276 include 'COMMON.SBRIDGE'
277 C Read bridging residues.
278 read (inp,*) ns,(iss(i),i=1,ns)
280 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
281 C Check whether the specified bridging residues are cystines.
283 if (itype(iss(i)).ne.1) then
284 write (iout,'(2a,i3,a)')
285 & 'Do you REALLY think that the residue ',
286 & restyp(itype(iss(i))),i,
287 & ' can form a disulfide bridge?!!!'
288 write (*,'(2a,i3,a)')
289 & 'Do you REALLY think that the residue ',
290 & restyp(itype(iss(i))),i,
291 & ' can form a disulfide bridge?!!!'
295 C Read preformed bridges.
297 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
298 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
301 C Check if the residues involved in bridges are in the specified list of
305 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
306 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
307 write (iout,'(a,i3,a)') 'Disulfide pair',i,
308 & ' contains residues present in other pairs.'
309 write (*,'(a,i3,a)') 'Disulfide pair',i,
310 & ' contains residues present in other pairs.'
315 if (ihpb(i).eq.iss(j)) goto 10
317 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
320 if (jhpb(i).eq.iss(j)) goto 20
322 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
335 c------------------------------------------------------------------------------
336 subroutine read_angles(kanal,iscor,energ,iprot,*)
337 implicit real*8 (a-h,o-z)
339 include 'DIMENSIONS.ZSCOPT'
340 include 'COMMON.INTERACT'
341 include 'COMMON.SBRIDGE'
344 include 'COMMON.CHAIN'
345 include 'COMMON.IOUNITS'
347 read(kanal,'(a80)',end=10,err=10) lineh
348 read(lineh(:5),*,err=8) ic
349 read(lineh(6:),*,err=8) energ
352 print *,'error, assuming e=1d10',lineh
356 read(lineh(18:),*,end=10,err=10) nss
358 read (lineh(20:),*,end=10,err=10)
359 & (IHPB(I),JHPB(I),I=1,NSS),iscor
361 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
362 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
365 c print *,"energy",energ," iscor",iscor
366 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
367 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
368 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
369 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
371 theta(i)=deg2rad*theta(i)
372 phi(i)=deg2rad*phi(i)
373 alph(i)=deg2rad*alph(i)
374 omeg(i)=deg2rad*omeg(i)
379 c-------------------------------------------------------------------------------
380 subroutine read_saxs_constr
381 implicit real*8 (a-h,o-z)
383 include 'DIMENSIONS.ZSCOPT'
384 include 'DIMENSIONS.FREE'
388 include 'COMMON.CONTROL'
389 include 'COMMON.CHAIN'
390 include 'COMMON.IOUNITS'
391 include 'COMMON.SBRIDGE'
392 double precision cm(3)
394 write (iout,*) "Calling read_saxs nsaxs",nsaxs
396 if (saxs_mode.eq.0) then
397 c SAXS distance distribution
399 read(inp,*) distsaxs(i),Psaxs(i)
403 Cnorm = Cnorm + Psaxs(i)
405 write (iout,*) "Cnorm",Cnorm
407 Psaxs(i)=Psaxs(i)/Cnorm
409 write (iout,*) "Normalized distance distribution from SAXS"
411 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
415 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
417 write (iout,*) "Wsaxs0",Wsaxs0
421 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
428 cm(j)=cm(j)+Csaxs(j,i)
436 Csaxs(j,i)=Csaxs(j,i)-cm(j)
439 write (iout,*) "SAXS sphere coordinates"
441 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)