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 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1233 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1234 if (homol_nset.gt.1)then
1235 call readi(controlcard,"ISET",iset,1)
1236 call card_concat(controlcard)
1237 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1240 waga_homology(1)=1.0
1244 write(iout,*) "read_constr_homology iset",iset
1245 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1252 cd write (iout,*) "nnt",nnt," nct",nct
1262 c Reading HM global scores (prob not required)
1265 do k=1,constr_homology
1269 c open (4,file="HMscore")
1270 c do k=1,constr_homology
1271 c read (4,*,end=521) hmscore_tmp
1272 c hmscore(k)=hmscore_tmp ! Another transformation can be used
1273 c write(*,*) "Model", k, ":", hmscore(k)
1284 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1286 if (read_homol_frag) then
1287 call read_klapaucjusz
1289 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1290 do k=1,constr_homology
1292 read(inp,'(a)') pdbfile
1293 c write (iout,*) "k ",k," pdbfile ",pdbfile
1294 c Next stament causes error upon compilation (?)
1295 c if(me.eq.king.or. .not. out1file)
1296 c write (iout,'(2a)') 'PDB data will be read from file ',
1297 c & pdbfile(:ilen(pdbfile))
1298 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1299 & pdbfile(:ilen(pdbfile))
1300 open(ipdbin,file=pdbfile,status='old',err=33)
1302 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1303 & pdbfile(:ilen(pdbfile))
1306 c print *,'Begin reading pdb data'
1308 c Files containing res sim or local scores (former containing sigmas)
1311 write(kic2,'(bz,i2.2)') k
1313 tpl_k_rescore="template"//kic2//".sco"
1316 call readpdb_template(k)
1319 crefjlee(j,i)=c(j,i)
1324 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1325 & (crefjlee(j,i+nres),j=1,3)
1327 write (iout,*) "READ HOMOLOGY INFO"
1328 write (iout,*) "read_constr_homology x: after reading pdb file"
1329 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1330 write (iout,*) "waga_dist",waga_dist
1331 write (iout,*) "waga_angle",waga_angle
1332 write (iout,*) "waga_theta",waga_theta
1333 write (iout,*) "waga_d",waga_d
1334 write (iout,*) "dist_cut",dist_cut
1343 c Distance restraints
1346 C Copy the coordinates from reference coordinates (?)
1350 c write (iout,*) "c(",j,i,") =",c(j,i)
1354 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1356 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1357 open (ientin,file=tpl_k_rescore,status='old')
1358 if (nnt.gt.1) rescore(k,1)=0.0d0
1359 do irec=nnt,nct ! loop for reading res sim
1360 if (read2sigma) then
1361 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1362 & rescore3_tmp,idomain_tmp
1364 idomain(k,i_tmp)=idomain_tmp
1365 rescore(k,i_tmp)=rescore_tmp
1366 rescore2(k,i_tmp)=rescore2_tmp
1367 rescore3(k,i_tmp)=rescore3_tmp
1368 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1369 & i_tmp,rescore2_tmp,rescore_tmp,
1370 & rescore3_tmp,idomain_tmp
1373 read (ientin,*,end=1401) rescore_tmp
1375 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1376 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1377 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1382 if (waga_dist.ne.0.0d0) then
1390 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1391 c write (iout,*) k,i,j,distal,dist2_cut
1393 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1394 & .and. distal.le.dist2_cut ) then
1400 c write (iout,*) "k",k
1401 c write (iout,*) "i",i," j",j," constr_homology",
1406 if (read2sigma) then
1409 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1411 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1412 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1413 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1415 if (odl(k,ii).le.dist_cut) then
1416 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1419 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1420 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1422 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1423 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1427 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1430 l_homo(k,ii)=.false.
1437 c Theta, dihedral and SC retraints
1439 if (waga_angle.gt.0.0d0) then
1440 c open (ientin,file=tpl_k_sigma_dih,status='old')
1441 c do irec=1,maxres-3 ! loop for reading sigma_dih
1442 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1443 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1444 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1445 c & sigma_dih(k,i+nnt-1)
1450 if (idomain(k,i).eq.0) then
1454 dih(k,i)=phiref(i) ! right?
1455 c read (ientin,*) sigma_dih(k,i) ! original variant
1456 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1457 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1458 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1459 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1460 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1462 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1463 & rescore(k,i-2)+rescore(k,i-3))/4.0
1464 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1465 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1466 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1467 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1468 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1469 c Instead of res sim other local measure of b/b str reliability possible
1470 if (sigma_dih(k,i).ne.0)
1471 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1472 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1477 if (waga_theta.gt.0.0d0) then
1478 c open (ientin,file=tpl_k_sigma_theta,status='old')
1479 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1480 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1481 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1482 c & sigma_theta(k,i+nnt-1)
1487 do i = nnt+2,nct ! right? without parallel.
1488 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1489 c do i=ithet_start,ithet_end ! with FG parallel.
1490 if (idomain(k,i).eq.0) then
1491 sigma_theta(k,i)=0.0
1494 thetatpl(k,i)=thetaref(i)
1495 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1496 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1497 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1498 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1499 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1500 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1501 & rescore(k,i-2))/3.0
1502 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1503 if (sigma_theta(k,i).ne.0)
1504 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1506 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1507 c rescore(k,i-2) ! right expression ?
1508 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1512 if (waga_d.gt.0.0d0) then
1513 c open (ientin,file=tpl_k_sigma_d,status='old')
1514 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1515 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1516 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1517 c & sigma_d(k,i+nnt-1)
1521 do i = nnt,nct ! right? without parallel.
1522 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1523 c do i=loc_start,loc_end ! with FG parallel.
1524 if (itype(i).eq.10) cycle
1525 if (idomain(k,i).eq.0 ) then
1532 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1533 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1534 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1535 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1536 sigma_d(k,i)=rescore3(k,i) ! right expression ?
1537 if (sigma_d(k,i).ne.0)
1538 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1540 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1541 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1542 c read (ientin,*) sigma_d(k,i) ! 1st variant
1547 c remove distance restraints not used in any model from the list
1548 c shift data in all arrays
1550 if (waga_dist.ne.0.0d0) then
1556 if (ii_in_use(ii).eq.0.and.liiflag) then
1560 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1561 & .not.liiflag.and.ii.eq.lim_odl) then
1562 if (ii.eq.lim_odl) then
1563 iishift=ii-iistart+1
1568 do ki=iistart,lim_odl-iishift
1569 ires_homo(ki)=ires_homo(ki+iishift)
1570 jres_homo(ki)=jres_homo(ki+iishift)
1571 ii_in_use(ki)=ii_in_use(ki+iishift)
1572 do k=1,constr_homology
1573 odl(k,ki)=odl(k,ki+iishift)
1574 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1575 l_homo(k,ki)=l_homo(k,ki+iishift)
1579 lim_odl=lim_odl-iishift
1585 endif ! .not. klapaucjusz
1587 if (constr_homology.gt.0) call homology_partition
1588 if (constr_homology.gt.0) call init_int_table
1589 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1590 cd & "lim_xx=",lim_xx
1591 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1592 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1596 if (.not.lprn) return
1597 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1598 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1599 write (iout,*) "Distance restraints from templates"
1601 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1602 & ii,ires_homo(ii),jres_homo(ii),
1603 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1604 & ki=1,constr_homology)
1606 write (iout,*) "Dihedral angle restraints from templates"
1608 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1609 & (rad2deg*dih(ki,i),
1610 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1612 write (iout,*) "Virtual-bond angle restraints from templates"
1614 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1615 & (rad2deg*thetatpl(ki,i),
1616 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1618 write (iout,*) "SC restraints from templates"
1620 write(iout,'(i5,100(4f8.2,4x))') i,
1621 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1622 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1625 c -----------------------------------------------------------------
1628 c----------------------------------------------------------------------
1629 subroutine read_klapaucjusz
1631 include 'DIMENSIONS'
1635 include 'COMMON.SETUP'
1636 include 'COMMON.CONTROL'
1637 include 'COMMON.CHAIN'
1638 include 'COMMON.IOUNITS'
1639 include 'COMMON.GEO'
1640 include 'COMMON.INTERACT'
1641 include 'COMMON.NAMES'
1642 include 'COMMON.HOMRESTR'
1643 character*256 fragfile
1644 integer ninclust(maxclust),inclust(max_template,maxclust),
1645 & nresclust(maxclust),iresclust(maxres,maxclust)
1648 character*24 model_ki_dist, model_ki_angle
1649 character*500 controlcard
1650 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1651 integer idomain(max_template,maxres)
1652 logical lprn /.true./
1655 logical unres_pdb,liiflag
1658 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1659 double precision, dimension (max_template,maxres) :: rescore
1660 double precision, dimension (max_template,maxres) :: rescore2
1661 character*24 tpl_k_rescore
1664 c For new homol impl
1666 include 'COMMON.VAR'
1668 double precision chomo(3,maxres2+2,max_template)
1669 call getenv("FRAGFILE",fragfile)
1670 open(ientin,file=fragfile,status="old",err=10)
1671 read(ientin,*) constr_homology,nclust
1677 do k=1,constr_homology
1678 read(ientin,'(a)') pdbfile
1679 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1680 & pdbfile(:ilen(pdbfile))
1681 open(ipdbin,file=pdbfile,status='old',err=33)
1683 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1684 & pdbfile(:ilen(pdbfile))
1701 read(ientin,*) ninclust(i),nresclust(i)
1702 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1703 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1706 c Loop over clusters
1709 do ll = 1,ninclust(l)
1717 idomain(k,iresclust(i,l)+1) = 1
1719 idomain(k,iresclust(i,l)) = 1
1723 c Distance restraints
1726 C Copy the coordinates from reference coordinates (?)
1730 c write (iout,*) "c(",j,i,") =",c(j,i)
1733 call int_from_cart(.true.,.false.)
1734 call sc_loc_geom(.false.)
1736 thetaref(i)=theta(i)
1739 if (waga_dist.ne.0.0d0) then
1747 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1748 c write (iout,*) k,i,j,distal,dist2_cut
1750 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1751 & .and. distal.le.dist2_cut ) then
1757 c write (iout,*) "k",k
1758 c write (iout,*) "i",i," j",j," constr_homology",
1763 if (read2sigma) then
1766 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1768 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1769 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1770 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1772 if (odl(k,ii).le.dist_cut) then
1773 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1776 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1777 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1779 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1780 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1784 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1787 c l_homo(k,ii)=.false.
1794 c Theta, dihedral and SC retraints
1796 if (waga_angle.gt.0.0d0) then
1798 if (idomain(k,i).eq.0) then
1799 c sigma_dih(k,i)=0.0
1803 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1804 & rescore(k,i-2)+rescore(k,i-3))/4.0
1805 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1806 c & " sigma_dihed",sigma_dih(k,i)
1807 if (sigma_dih(k,i).ne.0)
1808 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1813 if (waga_theta.gt.0.0d0) then
1815 if (idomain(k,i).eq.0) then
1816 c sigma_theta(k,i)=0.0
1819 thetatpl(k,i)=thetaref(i)
1820 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1821 & rescore(k,i-2))/3.0
1822 if (sigma_theta(k,i).ne.0)
1823 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1827 if (waga_d.gt.0.0d0) then
1829 if (itype(i).eq.10) cycle
1830 if (idomain(k,i).eq.0 ) then
1837 sigma_d(k,i)=rescore(k,i)
1838 if (sigma_d(k,i).ne.0)
1839 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1840 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1846 c remove distance restraints not used in any model from the list
1847 c shift data in all arrays
1849 if (waga_dist.ne.0.0d0) then
1855 if (ii_in_use(ii).eq.0.and.liiflag) then
1859 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1860 & .not.liiflag.and.ii.eq.lim_odl) then
1861 if (ii.eq.lim_odl) then
1862 iishift=ii-iistart+1
1867 do ki=iistart,lim_odl-iishift
1868 ires_homo(ki)=ires_homo(ki+iishift)
1869 jres_homo(ki)=jres_homo(ki+iishift)
1870 ii_in_use(ki)=ii_in_use(ki+iishift)
1871 do k=1,constr_homology
1872 odl(k,ki)=odl(k,ki+iishift)
1873 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1874 l_homo(k,ki)=l_homo(k,ki+iishift)
1878 lim_odl=lim_odl-iishift
1885 10 stop "Error infragment file"
1887 c----------------------------------------------------------------------