1 subroutine read_control
8 include 'COMMON.IOUNITS'
10 include 'COMMON.SBRIDGE'
11 include 'COMMON.CONTROL'
12 include 'COMMON.CLUSTER'
13 include 'COMMON.CHAIN'
14 include 'COMMON.HEADER'
15 include 'COMMON.FFIELD'
16 include 'COMMON.INTERACT'
17 include "COMMON.SPLITELE"
18 include 'COMMON.SHIELD'
20 include 'COMMON.LANGEVIN'
21 character*320 controlcard,ucase
25 integer i,i1,i2,it1,it2
27 read (INP,'(a80)') titel
28 call card_concat(controlcard)
30 call readi(controlcard,'NRES',nres,0)
31 call readi(controlcard,'RESCALE',rescale_mode,2)
32 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
33 write (iout,*) "DISTCHAINMAX",distchainmax
34 C Reading the dimensions of box in x,y,z coordinates
35 call reada(controlcard,'BOXX',boxxsize,100.0d0)
36 call reada(controlcard,'BOXY',boxysize,100.0d0)
37 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
38 c Cutoff range for interactions
39 call reada(controlcard,"R_CUT",r_cut,15.0d0)
40 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
41 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
42 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
43 if (lipthick.gt.0.0d0) then
44 bordliptop=(boxzsize+lipthick)/2.0
45 bordlipbot=bordliptop-lipthick
47 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
48 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
49 buflipbot=bordlipbot+lipbufthick
50 bufliptop=bordliptop-lipbufthick
51 if ((lipbufthick*2.0d0).gt.lipthick)
52 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
54 write(iout,*) "bordliptop=",bordliptop
55 write(iout,*) "bordlipbot=",bordlipbot
56 write(iout,*) "bufliptop=",bufliptop
57 write(iout,*) "buflipbot=",buflipbot
59 call readi(controlcard,'SHIELD',shield_mode,0)
60 write (iout,*) "SHIELD MODE",shield_mode
61 if (shield_mode.gt.0) then
63 C VSolvSphere the volume of solving sphere
65 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
66 C there will be no distinction between proline peptide group and normal peptide
67 C group in case of shielding parameters
68 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
69 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
70 write (iout,*) VSolvSphere,VSolvSphere_div
71 C long axis of side chain
73 long_r_sidechain(i)=vbldsc0(1,i)
74 short_r_sidechain(i)=sigma0(i)
78 call readi(controlcard,'PDBOUT',outpdb,0)
79 call readi(controlcard,'MOL2OUT',outmol2,0)
80 refstr=(index(controlcard,'REFSTR').gt.0)
81 write (iout,*) "REFSTR",refstr
82 pdbref=(index(controlcard,'PDBREF').gt.0)
83 iscode=index(controlcard,'ONE_LETTER')
84 tree=(index(controlcard,'MAKE_TREE').gt.0)
85 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
86 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
87 write (iout,*) "with_dihed_constr ",with_dihed_constr,
88 & " CONSTR_DIST",constr_dist
89 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
90 write (iout,*) "with_theta_constr ",with_theta_constr
92 min_var=(index(controlcard,'MINVAR').gt.0)
93 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
94 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
95 call readi(controlcard,'NCUT',ncut,0)
96 call readi(controlcard,'NCLUST',nclust,5)
97 call readi(controlcard,'SYM',symetr,1)
98 write (iout,*) 'sym', symetr
99 call readi(controlcard,'NSTART',nstart,0)
100 call readi(controlcard,'NEND',nend,0)
101 call reada(controlcard,'ECUT',ecut,10.0d0)
102 call reada(controlcard,'PROB',prob_limit,0.99d0)
103 write (iout,*) "Probability limit",prob_limit
104 lgrp=(index(controlcard,'LGRP').gt.0)
105 caonly=(index(controlcard,'CA_ONLY').gt.0)
106 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
108 & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
109 call readi(controlcard,'IOPT',iopt,2)
110 lside = index(controlcard,"SIDE").gt.0
111 efree = index(controlcard,"EFREE").gt.0
112 call readi(controlcard,'NTEMP',nT,1)
113 write (iout,*) "nT",nT
114 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
115 write (iout,*) "nT",nT
116 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
118 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
120 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
121 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
122 lprint_int=index(controlcard,"PRINT_INT") .gt.0
123 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
124 write (iout,*) "with_homology_constr ",with_dihed_constr,
125 & " CONSTR_HOMOLOGY",constr_homology
126 print_homology_restraints=
127 & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
128 print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
129 print_homology_models=
130 & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
131 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
132 call readi(controlcard,'NSAXS',nsaxs,0)
133 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
134 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
135 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
136 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
137 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
141 c--------------------------------------------------------------------------
144 C Read molecular data.
148 include 'COMMON.IOUNITS'
151 include 'COMMON.INTERACT'
152 include 'COMMON.LOCAL'
153 include 'COMMON.NAMES'
154 include 'COMMON.CHAIN'
155 include 'COMMON.FFIELD'
156 include 'COMMON.SBRIDGE'
157 include 'COMMON.HEADER'
158 include 'COMMON.CONTROL'
159 include 'COMMON.CONTACTS'
160 include 'COMMON.TIME1'
161 include 'COMMON.TORCNSTR'
162 include 'COMMON.SHIELD'
164 include 'COMMON.INFO'
166 character*4 sequence(maxres)
167 character*800 weightcard
169 double precision x(maxvar)
170 integer itype_pdb(maxres)
172 integer i,j,kkk,i1,i2,it1,it2
176 C Read weights of the subsequent energy terms.
177 call card_concat(weightcard)
178 write(iout,*) weightcard
179 call reada(weightcard,'WSC',wsc,1.0d0)
180 call reada(weightcard,'WLONG',wsc,wsc)
181 call reada(weightcard,'WSCP',wscp,1.0d0)
182 call reada(weightcard,'WELEC',welec,1.0D0)
183 call reada(weightcard,'WVDWPP',wvdwpp,welec)
184 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
185 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
186 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
187 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
188 call reada(weightcard,'WTURN3',wturn3,1.0D0)
189 call reada(weightcard,'WTURN4',wturn4,1.0D0)
190 call reada(weightcard,'WTURN6',wturn6,1.0D0)
191 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
192 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
193 call reada(weightcard,'WBOND',wbond,1.0D0)
194 call reada(weightcard,'WTOR',wtor,1.0D0)
195 call reada(weightcard,'WTORD',wtor_d,1.0D0)
196 call reada(weightcard,'WANG',wang,1.0D0)
197 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
198 call reada(weightcard,'WSAXS',wsaxs,0.0D0)
199 call reada(weightcard,'SCAL14',scal14,0.4D0)
200 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
201 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
202 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
203 if (index(weightcard,'SOFT').gt.0) ipot=6
204 call reada(weightcard,"D0CM",d0cm,3.78d0)
205 call reada(weightcard,"AKCM",akcm,15.1d0)
206 call reada(weightcard,"AKTH",akth,11.0d0)
207 call reada(weightcard,"AKCT",akct,12.0d0)
208 call reada(weightcard,"V1SS",v1ss,-1.08d0)
209 call reada(weightcard,"V2SS",v2ss,7.61d0)
210 call reada(weightcard,"V3SS",v3ss,13.7d0)
211 call reada(weightcard,"EBR",ebr,-5.50D0)
212 call reada(weightcard,'WSHIELD',wshield,1.0d0)
213 write(iout,*) 'WSHIELD',wshield
214 call reada(weightcard,'WLT',wliptran,0.0D0)
215 call reada(weightcard,"ATRISS",atriss,0.301D0)
216 call reada(weightcard,"BTRISS",btriss,0.021D0)
217 call reada(weightcard,"CTRISS",ctriss,1.001D0)
218 call reada(weightcard,"DTRISS",dtriss,1.001D0)
219 write (iout,*) "ATRISS=", atriss
220 write (iout,*) "BTRISS=", btriss
221 write (iout,*) "CTRISS=", ctriss
222 write (iout,*) "DTRISS=", dtriss
223 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
225 dyn_ss_mask(i)=.false.
229 dyn_ssbond_ij(i,j)=1.0d300
232 call reada(weightcard,"HT",Ht,0.0D0)
234 ss_depth=ebr/wsc-0.25*eps(1,1)
235 Ht=Ht/wsc-0.25*eps(1,1)
236 akcm=akcm*wstrain/wsc
237 akth=akth*wstrain/wsc
238 akct=akct*wstrain/wsc
239 v1ss=v1ss*wstrain/wsc
240 v2ss=v2ss*wstrain/wsc
241 v3ss=v3ss*wstrain/wsc
243 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
245 write (iout,'(/a)') "Disulfide bridge parameters:"
246 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
247 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
248 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
249 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
252 C 12/1/95 Added weight for the multi-body term WCORR
253 call reada(weightcard,'WCORRH',wcorr,1.0D0)
254 if (wcorr4.gt.0.0d0) wcorr=wcorr4
274 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
275 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
276 & wturn4,wturn6,wsccor
277 10 format (/'Energy-term weights (unscaled):'//
278 & 'WSCC= ',f10.6,' (SC-SC)'/
279 & 'WSCP= ',f10.6,' (SC-p)'/
280 & 'WELEC= ',f10.6,' (p-p electr)'/
281 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
282 & 'WBOND= ',f10.6,' (stretching)'/
283 & 'WANG= ',f10.6,' (bending)'/
284 & 'WSCLOC= ',f10.6,' (SC local)'/
285 & 'WTOR= ',f10.6,' (torsional)'/
286 & 'WTORD= ',f10.6,' (double torsional)'/
287 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
288 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
289 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
290 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
291 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
292 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
293 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
294 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
295 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
297 if (wcorr4.gt.0.0d0) then
298 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
299 & 'between contact pairs of peptide groups'
300 write (iout,'(2(a,f5.3/))')
301 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
302 & 'Range of quenching the correlation terms:',2*delt_corr
303 else if (wcorr.gt.0.0d0) then
304 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
305 & 'between contact pairs of peptide groups'
307 write (iout,'(a,f8.3)')
308 & 'Scaling factor of 1,4 SC-p interactions:',scal14
309 write (iout,'(a,f8.3)')
310 & 'General scaling factor of SC-p interactions:',scalscp
311 r0_corr=cutoff_corr-delt_corr
313 aad(i,1)=scalscp*aad(i,1)
314 aad(i,2)=scalscp*aad(i,2)
315 bad(i,1)=scalscp*bad(i,1)
316 bad(i,2)=scalscp*bad(i,2)
320 print *,'indpdb=',indpdb,' pdbref=',pdbref
322 C Read sequence if not taken from the pdb file.
323 if (iscode.gt.0) then
324 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
326 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
328 C Convert sequence to numeric code
330 itype(i)=rescode(i,sequence(i),iscode)
333 print '(20i4)',(itype(i),i=1,nres)
337 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
339 if (itype(i).eq.ntyp1) then
343 else if (iabs(itype(i+1)).ne.20) then
345 else if (iabs(itype(i)).ne.20) then
352 write (iout,*) "ITEL"
354 write (iout,*) i,itype(i),itel(i)
357 print *,'Call Read_Bridge.'
359 C this fragment reads diheadral constrains
360 if (with_dihed_constr) then
362 read (inp,*) ndih_constr
363 if (ndih_constr.gt.0) then
365 C write (iout,*) 'FTORS',ftors
366 C ftors is the force constant for torsional quartic constrains
367 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
370 & 'There are',ndih_constr,' constraints on phi angles.'
372 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
376 phi0(i)=deg2rad*phi0(i)
377 drange(i)=deg2rad*drange(i)
379 endif ! endif ndif_constr.gt.0
380 endif ! with_dihed_constr
381 if (with_theta_constr) then
382 C with_theta_constr is keyword allowing for occurance of theta constrains
383 read (inp,*) ntheta_constr
384 C ntheta_constr is the number of theta constrains
385 if (ntheta_constr.gt.0) then
387 read (inp,*) (itheta_constr(i),theta_constr0(i),
388 & theta_drange(i),for_thet_constr(i),
390 C the above code reads from 1 to ntheta_constr
391 C itheta_constr(i) residue i for which is theta_constr
392 C theta_constr0 the global minimum value
393 C theta_drange is range for which there is no energy penalty
394 C for_thet_constr is the force constant for quartic energy penalty
396 C if(me.eq.king.or..not.out1file)then
398 & 'There are',ntheta_constr,' constraints on phi angles.'
400 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
406 theta_constr0(i)=deg2rad*theta_constr0(i)
407 theta_drange(i)=deg2rad*theta_drange(i)
409 C if(me.eq.king.or..not.out1file)
410 C & write (iout,*) 'FTORS',ftors
411 C do i=1,ntheta_constr
412 C ii = itheta_constr(i)
413 C thetabound(1,ii) = phi0(i)-drange(i)
414 C thetabound(2,ii) = phi0(i)+drange(i)
416 endif ! ntheta_constr.gt.0
417 endif! with_theta_constr
421 print *,'NNT=',NNT,' NCT=',NCT
422 if (itype(1).eq.ntyp1) nnt=2
423 if (itype(nres).eq.ntyp1) nct=nct-1
424 if (nstart.lt.nnt) nstart=nnt
425 if (nend.gt.nct .or. nend.eq.0) nend=nct
426 write (iout,*) "nstart",nstart," nend",nend
427 write (iout,*) "calling read_saxs_consrtr",nsaxs
428 if (nsaxs.gt.0) call read_saxs_constr
430 if (constr_homology.gt.0) then
431 call read_constr_homology(print_homology_restraints)
435 c read(inp,'(a)') pdbfile
436 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
437 c open(ipdbin,file=pdbfile,status='old',err=33)
439 c 33 write (iout,'(a)') 'Error opening PDB file.'
442 c print *,'Begin reading pdb data'
444 c print *,'Finished reading pdb data'
445 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
447 c itype_pdb(i)=itype(i)
450 c write (iout,'(a,i3)') 'nsup=',nsup
452 c if (nsup.le.(nct-nnt+1)) then
453 c do i=0,nct-nnt+1-nsup
454 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
460 c & 'Error - sequences to be superposed do not match.'
463 c do i=0,nsup-(nct-nnt+1)
464 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
466 c nstart_sup=nstart_sup+i
472 c & 'Error - sequences to be superposed do not match.'
475 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
476 c & ' nstart_seq=',nstart_seq
480 write (iout,*) "molread: REFSTR",refstr
482 if (.not.pdbref) then
483 call read_angles(inp,*38)
485 38 write (iout,'(a)') 'Error reading reference structure.'
487 call mp_stopall(Error_Msg)
489 stop 'Error reading reference structure'
502 call contact(.true.,ncont_ref,icont_ref)
505 C write (iout,'(/a,i3,a)')
506 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
507 write (iout,'(20i4)') (iss(i),i=1,ns)
509 write(iout,*)"Running with dynamic disulfide-bond formation"
511 write (iout,'(/a/)') 'Pre-formed links are:'
517 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
518 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
524 if (ns.gt.0.and.dyn_ss) then
528 forcon(i-nss)=forcon(i)
535 dyn_ss_mask(iss(i))=.true.
538 c Read distance restraints
539 if (constr_dist.gt.0) then
540 call read_dist_constr
545 c-----------------------------------------------------------------------------
546 logical function seq_comp(itypea,itypeb,length)
548 integer length,itypea(length),itypeb(length)
551 if (itypea(i).ne.itypeb(i)) then
559 c-----------------------------------------------------------------------------
560 subroutine read_bridge
561 C Read information about disulfide bridges.
564 include 'COMMON.IOUNITS'
567 include 'COMMON.INTERACT'
568 include 'COMMON.LOCAL'
569 include 'COMMON.NAMES'
570 include 'COMMON.CHAIN'
571 include 'COMMON.FFIELD'
572 include 'COMMON.SBRIDGE'
573 include 'COMMON.HEADER'
574 include 'COMMON.CONTROL'
575 include 'COMMON.TIME1'
577 include 'COMMON.INFO'
580 C Read bridging residues.
581 read (inp,*) ns,(iss(i),i=1,ns)
583 C Check whether the specified bridging residues are cystines.
585 if (itype(iss(i)).ne.1) then
586 write (iout,'(2a,i3,a)')
587 & 'Do you REALLY think that the residue ',
588 & restyp(itype(iss(i))),i,
589 & ' can form a disulfide bridge?!!!'
590 write (*,'(2a,i3,a)')
591 & 'Do you REALLY think that the residue ',
592 & restyp(itype(iss(i))),i,
593 & ' can form a disulfide bridge?!!!'
595 call mp_stopall(error_msg)
601 C Read preformed bridges.
603 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
606 C Check if the residues involved in bridges are in the specified list of
610 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
611 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
612 write (iout,'(a,i3,a)') 'Disulfide pair',i,
613 & ' contains residues present in other pairs.'
614 write (*,'(a,i3,a)') 'Disulfide pair',i,
615 & ' contains residues present in other pairs.'
617 call mp_stopall(error_msg)
624 if (ihpb(i).eq.iss(j)) goto 10
626 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
629 if (jhpb(i).eq.iss(j)) goto 20
631 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
644 c----------------------------------------------------------------------------
645 subroutine read_angles(kanal,*)
650 include 'COMMON.CHAIN'
651 include 'COMMON.IOUNITS'
653 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
654 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
655 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
656 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
658 theta(i)=deg2rad*theta(i)
659 phi(i)=deg2rad*phi(i)
660 alph(i)=deg2rad*alph(i)
661 omeg(i)=deg2rad*omeg(i)
666 c----------------------------------------------------------------------------
667 subroutine reada(rekord,lancuch,wartosc,default)
669 character*(*) rekord,lancuch
670 double precision wartosc,default
673 iread=index(rekord,lancuch)
678 iread=iread+ilen(lancuch)+1
679 read (rekord(iread:),*) wartosc
682 c----------------------------------------------------------------------------
683 subroutine multreada(rekord,lancuch,tablica,dim,default)
686 double precision tablica(dim),default
687 character*(*) rekord,lancuch
693 iread=index(rekord,lancuch)
694 if (iread.eq.0) return
695 iread=iread+ilen(lancuch)+1
696 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
699 c----------------------------------------------------------------------------
700 subroutine readi(rekord,lancuch,wartosc,default)
702 character*(*) rekord,lancuch
703 integer wartosc,default
706 iread=index(rekord,lancuch)
711 iread=iread+ilen(lancuch)+1
712 read (rekord(iread:),*) wartosc
715 C----------------------------------------------------------------------
716 subroutine multreadi(rekord,lancuch,tablica,dim,default)
719 integer tablica(dim),default
720 character*(*) rekord,lancuch
727 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
728 if (iread.eq.0) return
729 iread=iread+ilen(lancuch)+1
730 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
734 c----------------------------------------------------------------------------
735 subroutine card_concat(card)
737 include 'COMMON.IOUNITS'
739 character*80 karta,ucase
741 read (inp,'(a)') karta
744 do while (karta(80:80).eq.'&')
745 card=card(:ilen(card)+1)//karta(:79)
746 read (inp,'(a)') karta
749 card=card(:ilen(card)+1)//karta
752 c----------------------------------------------------------------------------
761 include 'COMMON.IOUNITS'
762 include 'COMMON.CONTROL'
763 integer lenpre,lenpot,ilen
765 character*16 cformat,cprint
767 integer lenint,lenout
768 call getenv('INPUT',prefix)
769 call getenv('OUTPUT',prefout)
770 call getenv('INTIN',prefintin)
771 call getenv('COORD',cformat)
772 call getenv('PRINTCOOR',cprint)
773 call getenv('SCRATCHDIR',scratchdir)
776 if (index(ucase(cformat),'CX').gt.0) then
783 lenint=ilen(prefintin)
784 C Get the names and open the input files
785 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
787 write (liczba,'(bz,i3.3)') me
788 outname=prefout(:lenout)//'_clust.out_'//liczba
790 outname=prefout(:lenout)//'_clust.out'
793 intinname=prefintin(:lenint)//'.bx'
794 else if (from_cx) then
795 intinname=prefintin(:lenint)//'.cx'
797 intinname=prefintin(:lenint)//'.int'
799 rmsname=prefintin(:lenint)//'.rms'
800 open (jplot,file=prefout(:ilen(prefout))//'.tex',
802 open (jrms,file=rmsname,status='unknown')
803 open(iout,file=outname,status='unknown')
804 C Get parameter filenames and open the parameter files.
805 call getenv('BONDPAR',bondname)
806 open (ibond,file=bondname,status='old')
807 call getenv('THETPAR',thetname)
808 open (ithep,file=thetname,status='old')
809 call getenv('ROTPAR',rotname)
810 open (irotam,file=rotname,status='old')
811 call getenv('TORPAR',torname)
812 open (itorp,file=torname,status='old')
813 call getenv('TORDPAR',tordname)
814 open (itordp,file=tordname,status='old')
815 call getenv('FOURIER',fouriername)
816 open (ifourier,file=fouriername,status='old')
817 call getenv('ELEPAR',elename)
818 open (ielep,file=elename,status='old')
819 call getenv('SIDEPAR',sidename)
820 open (isidep,file=sidename,status='old')
821 call getenv('SIDEP',sidepname)
822 open (isidep1,file=sidepname,status="old")
823 call getenv('SCCORPAR',sccorname)
824 open (isccor,file=sccorname,status="old")
825 call getenv('LIPTRANPAR',liptranname)
826 open (iliptranpar,file=liptranname,status='old')
829 C 8/9/01 In the newest version SCp interaction constants are read from a file
830 C Use -DOLDSCP to use hard-coded constants instead.
832 call getenv('SCPPAR',scpname)
833 open (iscpp,file=scpname,status='old')
837 subroutine read_dist_constr
838 implicit real*8 (a-h,o-z)
843 include 'COMMON.CONTROL'
844 include 'COMMON.CHAIN'
845 include 'COMMON.IOUNITS'
846 include 'COMMON.SBRIDGE'
847 integer ifrag_(2,100),ipair_(2,100)
848 double precision wfrag_(100),wpair_(100)
849 character*500 controlcard
850 logical lprn /.true./
851 logical normalize,next
853 double precision xlink(4,0:4) /
855 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
856 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
857 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
858 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
859 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
860 write (iout,*) "Calling read_dist_constr"
861 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
867 call card_concat(controlcard)
868 next = index(controlcard,"NEXT").gt.0
869 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
870 write (iout,*) "restr_type",restr_type
871 call readi(controlcard,"NFRAG",nfrag_,0)
872 call readi(controlcard,"NFRAG",nfrag_,0)
873 call readi(controlcard,"NPAIR",npair_,0)
874 call readi(controlcard,"NDIST",ndist_,0)
875 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
876 if (restr_type.eq.10)
877 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
878 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
879 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
880 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
881 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
882 normalize = index(controlcard,"NORMALIZE").gt.0
883 write (iout,*) "WBOLTZD",wboltzd
884 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
885 c write (iout,*) "IFRAG"
887 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
889 c write (iout,*) "IPAIR"
891 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
893 if (nfrag_.gt.0) then
895 read(inp,'(a)') pdbfile
897 & "Distance restraints will be constructed from structure ",pdbfile
898 open(ipdbin,file=pdbfile,status='old',err=11)
904 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
905 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
906 & ifrag_(2,i)=nstart_sup+nsup-1
907 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
909 if (wfrag_(i).eq.0.0d0) cycle
910 do j=ifrag_(1,i),ifrag_(2,i)-1
912 c write (iout,*) "j",j," k",k
914 if (restr_type.eq.1) then
920 forcon(nhpb)=wfrag_(i)
921 else if (constr_dist.eq.2) then
922 if (ddjk.le.dist_cut) then
928 forcon(nhpb)=wfrag_(i)
930 else if (restr_type.eq.3) then
936 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
938 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
939 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
944 if (wpair_(i).eq.0.0d0) cycle
952 do j=ifrag_(1,ii),ifrag_(2,ii)
953 do k=ifrag_(1,jj),ifrag_(2,jj)
954 if (restr_type.eq.1) then
960 forcon(nhpb)=wfrag_(i)
961 else if (constr_dist.eq.2) then
962 if (ddjk.le.dist_cut) then
968 forcon(nhpb)=wfrag_(i)
970 else if (restr_type.eq.3) then
976 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
978 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
979 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
985 write (iout,*) "Distance restraints as read from input"
987 if (restr_type.eq.11) then
988 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
989 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
990 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
991 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
994 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
995 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
996 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
997 if (ibecarb(nhpb).gt.0) then
998 ihpb(nhpb)=ihpb(nhpb)+nres
999 jhpb(nhpb)=jhpb(nhpb)+nres
1001 else if (constr_dist.eq.10) then
1002 c Cross-lonk Markov-like potential
1003 call card_concat(controlcard)
1004 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1005 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1007 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1008 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1009 if (index(controlcard,"ZL").gt.0) then
1011 else if (index(controlcard,"ADH").gt.0) then
1013 else if (index(controlcard,"PDH").gt.0) then
1015 else if (index(controlcard,"DSS").gt.0) then
1020 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1021 & xlink(1,link_type))
1022 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1023 & xlink(2,link_type))
1024 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1025 & xlink(3,link_type))
1026 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1027 & xlink(4,link_type))
1028 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1029 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1030 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1031 if (forcon(nhpb+1).le.0.0d0 .or.
1032 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1034 irestr_type(nhpb)=10
1035 if (ibecarb(nhpb).gt.0) then
1036 ihpb(nhpb)=ihpb(nhpb)+nres
1037 jhpb(nhpb)=jhpb(nhpb)+nres
1039 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1040 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1041 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1045 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1046 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1047 if (forcon(nhpb+1).gt.0.0d0) then
1049 if (dhpb1(nhpb).eq.0.0d0) then
1054 if (ibecarb(nhpb).gt.0) then
1055 ihpb(nhpb)=ihpb(nhpb)+nres
1056 jhpb(nhpb)=jhpb(nhpb)+nres
1058 if (dhpb(nhpb).eq.0.0d0)
1059 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1061 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1062 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1064 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1065 C if (forcon(nhpb+1).gt.0.0d0) then
1067 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1075 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1076 & fordepthmax=fordepth(i)
1079 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1082 if (nhpb.gt.nss) then
1083 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1084 & "The following",nhpb-nss,
1085 & " distance restraints have been imposed:",
1086 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1089 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1090 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1097 11 write (iout,*)"read_dist_restr: error reading reference structure"
1100 c-------------------------------------------------------------------------------
1101 subroutine read_saxs_constr
1102 implicit real*8 (a-h,o-z)
1103 include 'DIMENSIONS'
1107 include 'COMMON.CONTROL'
1108 include 'COMMON.CHAIN'
1109 include 'COMMON.IOUNITS'
1110 include 'COMMON.SBRIDGE'
1111 double precision cm(3)
1113 write (iout,*) "Calling read_saxs nsaxs",nsaxs
1115 if (saxs_mode.eq.0) then
1116 c SAXS distance distribution
1118 read(inp,*) distsaxs(i),Psaxs(i)
1122 Cnorm = Cnorm + Psaxs(i)
1124 write (iout,*) "Cnorm",Cnorm
1126 Psaxs(i)=Psaxs(i)/Cnorm
1128 write (iout,*) "Normalized distance distribution from SAXS"
1130 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1134 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1136 write (iout,*) "Wsaxs0",Wsaxs0
1140 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1147 cm(j)=cm(j)+Csaxs(j,i)
1155 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1158 write (iout,*) "SAXS sphere coordinates"
1160 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
1165 c====-------------------------------------------------------------------
1166 subroutine read_constr_homology(lprn)
1168 include 'DIMENSIONS'
1172 include 'COMMON.SETUP'
1173 include 'COMMON.CONTROL'
1174 include 'COMMON.CHAIN'
1175 include 'COMMON.IOUNITS'
1176 include 'COMMON.GEO'
1177 include 'COMMON.INTERACT'
1178 include 'COMMON.NAMES'
1179 include 'COMMON.HOMRESTR'
1181 c For new homol impl
1183 include 'COMMON.VAR'
1184 c include 'include_unres/COMMON.VAR'
1187 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1189 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1190 c & sigma_odl_temp(maxres,maxres,max_template)
1192 character*24 model_ki_dist, model_ki_angle
1193 character*500 controlcard
1194 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1195 integer idomain(max_template,maxres)
1199 logical unres_pdb,liiflag
1201 c FP - Nov. 2014 Temporary specifications for new vars
1203 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
1205 double precision, dimension (max_template,maxres) :: rescore
1206 double precision, dimension (max_template,maxres) :: rescore2
1207 double precision, dimension (max_template,maxres) :: rescore3
1208 character*24 tpl_k_rescore
1209 c -----------------------------------------------------------------
1210 c Reading multiple PDB ref structures and calculation of retraints
1211 c not using pre-computed ones stored in files model_ki_{dist,angle}
1213 c -----------------------------------------------------------------
1216 c Alternative: reading from input
1218 write (iout,*) "BEGIN READ HOMOLOGY INFO"
1225 call card_concat(controlcard)
1226 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1227 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1228 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1229 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1230 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1231 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1232 dist1cut=(index(controlcard,'DIST1CUT').gt.0)
1233 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1234 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1235 if (homol_nset.gt.1)then
1236 call readi(controlcard,"ISET",iset,1)
1237 call card_concat(controlcard)
1238 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1241 waga_homology(1)=1.0
1245 write(iout,*) "read_constr_homology iset",iset
1246 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1253 cd write (iout,*) "nnt",nnt," nct",nct
1263 c Reading HM global scores (prob not required)
1266 do k=1,constr_homology
1270 c open (4,file="HMscore")
1271 c do k=1,constr_homology
1272 c read (4,*,end=521) hmscore_tmp
1273 c hmscore(k)=hmscore_tmp ! Another transformation can be used
1274 c write(*,*) "Model", k, ":", hmscore(k)
1285 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1287 if (read_homol_frag) then
1288 call read_klapaucjusz
1290 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1291 do k=1,constr_homology
1293 read(inp,'(a)') pdbfile
1294 c write (iout,*) "k ",k," pdbfile ",pdbfile
1295 c Next stament causes error upon compilation (?)
1296 c if(me.eq.king.or. .not. out1file)
1297 c write (iout,'(2a)') 'PDB data will be read from file ',
1298 c & pdbfile(:ilen(pdbfile))
1299 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1300 & pdbfile(:ilen(pdbfile))
1301 open(ipdbin,file=pdbfile,status='old',err=33)
1303 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1304 & pdbfile(:ilen(pdbfile))
1307 c print *,'Begin reading pdb data'
1309 c Files containing res sim or local scores (former containing sigmas)
1312 write(kic2,'(bz,i2.2)') k
1314 tpl_k_rescore="template"//kic2//".sco"
1317 call readpdb_template(k)
1320 crefjlee(j,i)=c(j,i)
1325 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1326 & (crefjlee(j,i+nres),j=1,3)
1328 write (iout,*) "READ HOMOLOGY INFO"
1329 write (iout,*) "read_constr_homology x: after reading pdb file"
1330 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1331 write (iout,*) "waga_dist",waga_dist
1332 write (iout,*) "waga_angle",waga_angle
1333 write (iout,*) "waga_theta",waga_theta
1334 write (iout,*) "waga_d",waga_d
1335 write (iout,*) "dist_cut",dist_cut
1344 c Distance restraints
1347 C Copy the coordinates from reference coordinates (?)
1351 c write (iout,*) "c(",j,i,") =",c(j,i)
1355 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1357 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1358 open (ientin,file=tpl_k_rescore,status='old')
1359 if (nnt.gt.1) rescore(k,1)=0.0d0
1360 do irec=nnt,nct ! loop for reading res sim
1361 if (read2sigma) then
1362 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1363 & rescore3_tmp,idomain_tmp
1365 idomain(k,i_tmp)=idomain_tmp
1366 rescore(k,i_tmp)=rescore_tmp
1367 rescore2(k,i_tmp)=rescore2_tmp
1368 rescore3(k,i_tmp)=rescore3_tmp
1369 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1370 & i_tmp,rescore2_tmp,rescore_tmp,
1371 & rescore3_tmp,idomain_tmp
1374 read (ientin,*,end=1401) rescore_tmp
1376 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1377 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1378 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1383 if (waga_dist.ne.0.0d0) then
1391 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1392 c write (iout,*) k,i,j,distal,dist2_cut
1393 if (dist1cut .and. k.gt.1) then
1395 if (l_homo(1,ii)) then
1401 sigma_odl(k,ii)=sigma_odl(1,ii)
1403 l_homo(k,ii)=.false.
1406 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1407 & .and. distal.le.dist2_cut ) then
1413 c write (iout,*) "k",k
1414 c write (iout,*) "i",i," j",j," constr_homology",
1419 if (read2sigma) then
1422 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1424 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1425 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1426 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1428 if (odl(k,ii).le.dist_cut) then
1429 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1432 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1433 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1435 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1436 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1440 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1443 l_homo(k,ii)=.false.
1451 c Theta, dihedral and SC retraints
1453 if (waga_angle.gt.0.0d0) then
1454 c open (ientin,file=tpl_k_sigma_dih,status='old')
1455 c do irec=1,maxres-3 ! loop for reading sigma_dih
1456 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1457 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1458 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1459 c & sigma_dih(k,i+nnt-1)
1464 if (idomain(k,i).eq.0) then
1468 dih(k,i)=phiref(i) ! right?
1469 c read (ientin,*) sigma_dih(k,i) ! original variant
1470 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1471 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1472 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1473 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1474 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1476 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1477 & rescore(k,i-2)+rescore(k,i-3))/4.0
1478 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1479 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1480 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1481 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1482 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1483 c Instead of res sim other local measure of b/b str reliability possible
1484 if (sigma_dih(k,i).ne.0)
1485 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1486 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1491 if (waga_theta.gt.0.0d0) then
1492 c open (ientin,file=tpl_k_sigma_theta,status='old')
1493 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1494 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1495 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1496 c & sigma_theta(k,i+nnt-1)
1501 do i = nnt+2,nct ! right? without parallel.
1502 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1503 c do i=ithet_start,ithet_end ! with FG parallel.
1504 if (idomain(k,i).eq.0) then
1505 sigma_theta(k,i)=0.0
1508 thetatpl(k,i)=thetaref(i)
1509 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1510 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1511 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1512 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1513 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1514 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1515 & rescore(k,i-2))/3.0
1516 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1517 if (sigma_theta(k,i).ne.0)
1518 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1520 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1521 c rescore(k,i-2) ! right expression ?
1522 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1526 if (waga_d.gt.0.0d0) then
1527 c open (ientin,file=tpl_k_sigma_d,status='old')
1528 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1529 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1530 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1531 c & sigma_d(k,i+nnt-1)
1535 do i = nnt,nct ! right? without parallel.
1536 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1537 c do i=loc_start,loc_end ! with FG parallel.
1538 if (itype(i).eq.10) cycle
1539 if (idomain(k,i).eq.0 ) then
1546 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1547 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1548 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1549 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1550 sigma_d(k,i)=rescore3(k,i) ! right expression ?
1551 if (sigma_d(k,i).ne.0)
1552 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1554 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1555 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1556 c read (ientin,*) sigma_d(k,i) ! 1st variant
1561 c remove distance restraints not used in any model from the list
1562 c shift data in all arrays
1564 if (waga_dist.ne.0.0d0) then
1570 if (ii_in_use(ii).eq.0.and.liiflag) then
1574 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1575 & .not.liiflag.and.ii.eq.lim_odl) then
1576 if (ii.eq.lim_odl) then
1577 iishift=ii-iistart+1
1582 do ki=iistart,lim_odl-iishift
1583 ires_homo(ki)=ires_homo(ki+iishift)
1584 jres_homo(ki)=jres_homo(ki+iishift)
1585 ii_in_use(ki)=ii_in_use(ki+iishift)
1586 do k=1,constr_homology
1587 odl(k,ki)=odl(k,ki+iishift)
1588 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1589 l_homo(k,ki)=l_homo(k,ki+iishift)
1593 lim_odl=lim_odl-iishift
1599 endif ! .not. klapaucjusz
1601 if (constr_homology.gt.0) call homology_partition
1602 if (constr_homology.gt.0) call init_int_table
1603 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1604 cd & "lim_xx=",lim_xx
1605 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1606 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1610 if (.not.lprn) return
1611 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1612 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1613 write (iout,*) "Distance restraints from templates"
1615 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1616 & ii,ires_homo(ii),jres_homo(ii),
1617 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1618 & ki=1,constr_homology)
1620 write (iout,*) "Dihedral angle restraints from templates"
1622 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1623 & (rad2deg*dih(ki,i),
1624 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1626 write (iout,*) "Virtual-bond angle restraints from templates"
1628 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1629 & (rad2deg*thetatpl(ki,i),
1630 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1632 write (iout,*) "SC restraints from templates"
1634 write(iout,'(i5,100(4f8.2,4x))') i,
1635 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1636 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1639 c -----------------------------------------------------------------
1642 c----------------------------------------------------------------------
1643 subroutine read_klapaucjusz
1645 include 'DIMENSIONS'
1649 include 'COMMON.SETUP'
1650 include 'COMMON.CONTROL'
1651 include 'COMMON.CHAIN'
1652 include 'COMMON.IOUNITS'
1653 include 'COMMON.GEO'
1654 include 'COMMON.INTERACT'
1655 include 'COMMON.NAMES'
1656 include 'COMMON.HOMRESTR'
1657 character*256 fragfile
1658 integer ninclust(maxclust),inclust(max_template,maxclust),
1659 & nresclust(maxclust),iresclust(maxres,maxclust)
1662 character*24 model_ki_dist, model_ki_angle
1663 character*500 controlcard
1664 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1665 integer idomain(max_template,maxres)
1666 logical lprn /.true./
1669 logical unres_pdb,liiflag
1672 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1673 double precision, dimension (max_template,maxres) :: rescore
1674 double precision, dimension (max_template,maxres) :: rescore2
1675 character*24 tpl_k_rescore
1678 c For new homol impl
1680 include 'COMMON.VAR'
1682 double precision chomo(3,maxres2+2,max_template)
1683 call getenv("FRAGFILE",fragfile)
1684 open(ientin,file=fragfile,status="old",err=10)
1685 read(ientin,*) constr_homology,nclust
1691 do k=1,constr_homology
1692 read(ientin,'(a)') pdbfile
1693 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1694 & pdbfile(:ilen(pdbfile))
1695 open(ipdbin,file=pdbfile,status='old',err=33)
1697 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1698 & pdbfile(:ilen(pdbfile))
1702 call readpdb_template(k)
1715 read(ientin,*) ninclust(i),nresclust(i)
1716 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1717 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1720 c Loop over clusters
1723 do ll = 1,ninclust(l)
1731 idomain(k,iresclust(i,l)+1) = 1
1733 idomain(k,iresclust(i,l)) = 1
1737 c Distance restraints
1740 C Copy the coordinates from reference coordinates (?)
1744 c write (iout,*) "c(",j,i,") =",c(j,i)
1747 call int_from_cart1(.false.)
1748 call int_from_cart(.true.,.false.)
1749 call sc_loc_geom(.false.)
1751 thetaref(i)=theta(i)
1754 if (waga_dist.ne.0.0d0) then
1762 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1763 c write (iout,*) k,i,j,distal,dist2_cut
1765 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1766 & .and. distal.le.dist2_cut ) then
1772 c write (iout,*) "k",k
1773 c write (iout,*) "i",i," j",j," constr_homology",
1778 if (read2sigma) then
1781 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1783 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1784 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1785 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1787 if (odl(k,ii).le.dist_cut) then
1788 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1791 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1792 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1794 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1795 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1799 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1802 c l_homo(k,ii)=.false.
1809 c Theta, dihedral and SC retraints
1811 if (waga_angle.gt.0.0d0) then
1813 if (idomain(k,i).eq.0) then
1814 c sigma_dih(k,i)=0.0
1818 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1819 & rescore(k,i-2)+rescore(k,i-3))/4.0
1820 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1821 c & " sigma_dihed",sigma_dih(k,i)
1822 if (sigma_dih(k,i).ne.0)
1823 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1828 if (waga_theta.gt.0.0d0) then
1830 if (idomain(k,i).eq.0) then
1831 c sigma_theta(k,i)=0.0
1834 thetatpl(k,i)=thetaref(i)
1835 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1836 & rescore(k,i-2))/3.0
1837 if (sigma_theta(k,i).ne.0)
1838 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1842 if (waga_d.gt.0.0d0) then
1844 if (itype(i).eq.10) cycle
1845 if (idomain(k,i).eq.0 ) then
1852 sigma_d(k,i)=rescore(k,i)
1853 if (sigma_d(k,i).ne.0)
1854 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1855 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1861 c remove distance restraints not used in any model from the list
1862 c shift data in all arrays
1864 if (waga_dist.ne.0.0d0) then
1870 if (ii_in_use(ii).eq.0.and.liiflag) then
1874 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1875 & .not.liiflag.and.ii.eq.lim_odl) then
1876 if (ii.eq.lim_odl) then
1877 iishift=ii-iistart+1
1882 do ki=iistart,lim_odl-iishift
1883 ires_homo(ki)=ires_homo(ki+iishift)
1884 jres_homo(ki)=jres_homo(ki+iishift)
1885 ii_in_use(ki)=ii_in_use(ki+iishift)
1886 do k=1,constr_homology
1887 odl(k,ki)=odl(k,ki+iishift)
1888 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1889 l_homo(k,ki)=l_homo(k,ki+iishift)
1893 lim_odl=lim_odl-iishift
1900 10 stop "Error infragment file"
1902 c----------------------------------------------------------------------