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 call readi(controlcard,'NSAXS',nsaxs,0)
132 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
133 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
134 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
135 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
136 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
140 c--------------------------------------------------------------------------
143 C Read molecular data.
147 include 'COMMON.IOUNITS'
150 include 'COMMON.INTERACT'
151 include 'COMMON.LOCAL'
152 include 'COMMON.NAMES'
153 include 'COMMON.CHAIN'
154 include 'COMMON.FFIELD'
155 include 'COMMON.SBRIDGE'
156 include 'COMMON.HEADER'
157 include 'COMMON.CONTROL'
158 include 'COMMON.CONTACTS'
159 include 'COMMON.TIME1'
160 include 'COMMON.TORCNSTR'
161 include 'COMMON.SHIELD'
163 include 'COMMON.INFO'
165 character*4 sequence(maxres)
166 character*800 weightcard
168 double precision x(maxvar)
169 integer itype_pdb(maxres)
171 integer i,j,kkk,i1,i2,it1,it2
175 C Read weights of the subsequent energy terms.
176 call card_concat(weightcard)
177 write(iout,*) weightcard
178 call reada(weightcard,'WSC',wsc,1.0d0)
179 call reada(weightcard,'WLONG',wsc,wsc)
180 call reada(weightcard,'WSCP',wscp,1.0d0)
181 call reada(weightcard,'WELEC',welec,1.0D0)
182 call reada(weightcard,'WVDWPP',wvdwpp,welec)
183 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
184 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
185 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
186 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
187 call reada(weightcard,'WTURN3',wturn3,1.0D0)
188 call reada(weightcard,'WTURN4',wturn4,1.0D0)
189 call reada(weightcard,'WTURN6',wturn6,1.0D0)
190 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
191 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
192 call reada(weightcard,'WBOND',wbond,1.0D0)
193 call reada(weightcard,'WTOR',wtor,1.0D0)
194 call reada(weightcard,'WTORD',wtor_d,1.0D0)
195 call reada(weightcard,'WANG',wang,1.0D0)
196 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
197 call reada(weightcard,'WSAXS',wsaxs,0.0D0)
198 call reada(weightcard,'SCAL14',scal14,0.4D0)
199 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
200 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
201 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
202 if (index(weightcard,'SOFT').gt.0) ipot=6
203 call reada(weightcard,"D0CM",d0cm,3.78d0)
204 call reada(weightcard,"AKCM",akcm,15.1d0)
205 call reada(weightcard,"AKTH",akth,11.0d0)
206 call reada(weightcard,"AKCT",akct,12.0d0)
207 call reada(weightcard,"V1SS",v1ss,-1.08d0)
208 call reada(weightcard,"V2SS",v2ss,7.61d0)
209 call reada(weightcard,"V3SS",v3ss,13.7d0)
210 call reada(weightcard,"EBR",ebr,-5.50D0)
211 call reada(weightcard,'WSHIELD',wshield,1.0d0)
212 write(iout,*) 'WSHIELD',wshield
213 call reada(weightcard,'WLT',wliptran,0.0D0)
214 call reada(weightcard,"ATRISS",atriss,0.301D0)
215 call reada(weightcard,"BTRISS",btriss,0.021D0)
216 call reada(weightcard,"CTRISS",ctriss,1.001D0)
217 call reada(weightcard,"DTRISS",dtriss,1.001D0)
218 write (iout,*) "ATRISS=", atriss
219 write (iout,*) "BTRISS=", btriss
220 write (iout,*) "CTRISS=", ctriss
221 write (iout,*) "DTRISS=", dtriss
222 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
224 dyn_ss_mask(i)=.false.
228 dyn_ssbond_ij(i,j)=1.0d300
231 call reada(weightcard,"HT",Ht,0.0D0)
233 ss_depth=ebr/wsc-0.25*eps(1,1)
234 Ht=Ht/wsc-0.25*eps(1,1)
235 akcm=akcm*wstrain/wsc
236 akth=akth*wstrain/wsc
237 akct=akct*wstrain/wsc
238 v1ss=v1ss*wstrain/wsc
239 v2ss=v2ss*wstrain/wsc
240 v3ss=v3ss*wstrain/wsc
242 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
244 write (iout,'(/a)') "Disulfide bridge parameters:"
245 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
246 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
247 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
248 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
251 C 12/1/95 Added weight for the multi-body term WCORR
252 call reada(weightcard,'WCORRH',wcorr,1.0D0)
253 if (wcorr4.gt.0.0d0) wcorr=wcorr4
273 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
274 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
275 & wturn4,wturn6,wsccor
276 10 format (/'Energy-term weights (unscaled):'//
277 & 'WSCC= ',f10.6,' (SC-SC)'/
278 & 'WSCP= ',f10.6,' (SC-p)'/
279 & 'WELEC= ',f10.6,' (p-p electr)'/
280 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
281 & 'WBOND= ',f10.6,' (stretching)'/
282 & 'WANG= ',f10.6,' (bending)'/
283 & 'WSCLOC= ',f10.6,' (SC local)'/
284 & 'WTOR= ',f10.6,' (torsional)'/
285 & 'WTORD= ',f10.6,' (double torsional)'/
286 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
287 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
288 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
289 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
290 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
291 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
292 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
293 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
294 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
296 if (wcorr4.gt.0.0d0) then
297 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
298 & 'between contact pairs of peptide groups'
299 write (iout,'(2(a,f5.3/))')
300 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
301 & 'Range of quenching the correlation terms:',2*delt_corr
302 else if (wcorr.gt.0.0d0) then
303 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
304 & 'between contact pairs of peptide groups'
306 write (iout,'(a,f8.3)')
307 & 'Scaling factor of 1,4 SC-p interactions:',scal14
308 write (iout,'(a,f8.3)')
309 & 'General scaling factor of SC-p interactions:',scalscp
310 r0_corr=cutoff_corr-delt_corr
312 aad(i,1)=scalscp*aad(i,1)
313 aad(i,2)=scalscp*aad(i,2)
314 bad(i,1)=scalscp*bad(i,1)
315 bad(i,2)=scalscp*bad(i,2)
319 print *,'indpdb=',indpdb,' pdbref=',pdbref
321 C Read sequence if not taken from the pdb file.
322 if (iscode.gt.0) then
323 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
325 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
327 C Convert sequence to numeric code
329 itype(i)=rescode(i,sequence(i),iscode)
332 print '(20i4)',(itype(i),i=1,nres)
336 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
338 if (itype(i).eq.ntyp1) then
342 else if (iabs(itype(i+1)).ne.20) then
344 else if (iabs(itype(i)).ne.20) then
351 write (iout,*) "ITEL"
353 write (iout,*) i,itype(i),itel(i)
356 print *,'Call Read_Bridge.'
358 C this fragment reads diheadral constrains
359 if (with_dihed_constr) then
361 read (inp,*) ndih_constr
362 if (ndih_constr.gt.0) then
364 C write (iout,*) 'FTORS',ftors
365 C ftors is the force constant for torsional quartic constrains
366 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
369 & 'There are',ndih_constr,' constraints on phi angles.'
371 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
375 phi0(i)=deg2rad*phi0(i)
376 drange(i)=deg2rad*drange(i)
378 endif ! endif ndif_constr.gt.0
379 endif ! with_dihed_constr
380 if (with_theta_constr) then
381 C with_theta_constr is keyword allowing for occurance of theta constrains
382 read (inp,*) ntheta_constr
383 C ntheta_constr is the number of theta constrains
384 if (ntheta_constr.gt.0) then
386 read (inp,*) (itheta_constr(i),theta_constr0(i),
387 & theta_drange(i),for_thet_constr(i),
389 C the above code reads from 1 to ntheta_constr
390 C itheta_constr(i) residue i for which is theta_constr
391 C theta_constr0 the global minimum value
392 C theta_drange is range for which there is no energy penalty
393 C for_thet_constr is the force constant for quartic energy penalty
395 C if(me.eq.king.or..not.out1file)then
397 & 'There are',ntheta_constr,' constraints on phi angles.'
399 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
405 theta_constr0(i)=deg2rad*theta_constr0(i)
406 theta_drange(i)=deg2rad*theta_drange(i)
408 C if(me.eq.king.or..not.out1file)
409 C & write (iout,*) 'FTORS',ftors
410 C do i=1,ntheta_constr
411 C ii = itheta_constr(i)
412 C thetabound(1,ii) = phi0(i)-drange(i)
413 C thetabound(2,ii) = phi0(i)+drange(i)
415 endif ! ntheta_constr.gt.0
416 endif! with_theta_constr
420 print *,'NNT=',NNT,' NCT=',NCT
421 if (itype(1).eq.ntyp1) nnt=2
422 if (itype(nres).eq.ntyp1) nct=nct-1
423 if (nstart.lt.nnt) nstart=nnt
424 if (nend.gt.nct .or. nend.eq.0) nend=nct
425 write (iout,*) "nstart",nstart," nend",nend
426 write (iout,*) "calling read_saxs_consrtr",nsaxs
427 if (nsaxs.gt.0) call read_saxs_constr
429 if (constr_homology.gt.0) then
430 call read_constr_homology(print_homology_restraints)
434 c read(inp,'(a)') pdbfile
435 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
436 c open(ipdbin,file=pdbfile,status='old',err=33)
438 c 33 write (iout,'(a)') 'Error opening PDB file.'
441 c print *,'Begin reading pdb data'
443 c print *,'Finished reading pdb data'
444 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
446 c itype_pdb(i)=itype(i)
449 c write (iout,'(a,i3)') 'nsup=',nsup
451 c if (nsup.le.(nct-nnt+1)) then
452 c do i=0,nct-nnt+1-nsup
453 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
459 c & 'Error - sequences to be superposed do not match.'
462 c do i=0,nsup-(nct-nnt+1)
463 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
465 c nstart_sup=nstart_sup+i
471 c & 'Error - sequences to be superposed do not match.'
474 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
475 c & ' nstart_seq=',nstart_seq
479 write (iout,*) "molread: REFSTR",refstr
481 if (.not.pdbref) then
482 call read_angles(inp,*38)
484 38 write (iout,'(a)') 'Error reading reference structure.'
486 call mp_stopall(Error_Msg)
488 stop 'Error reading reference structure'
501 call contact(.true.,ncont_ref,icont_ref)
504 C write (iout,'(/a,i3,a)')
505 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
506 write (iout,'(20i4)') (iss(i),i=1,ns)
508 write(iout,*)"Running with dynamic disulfide-bond formation"
510 write (iout,'(/a/)') 'Pre-formed links are:'
516 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
517 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
523 if (ns.gt.0.and.dyn_ss) then
527 forcon(i-nss)=forcon(i)
534 dyn_ss_mask(iss(i))=.true.
537 c Read distance restraints
538 if (constr_dist.gt.0) then
539 call read_dist_constr
544 c-----------------------------------------------------------------------------
545 logical function seq_comp(itypea,itypeb,length)
547 integer length,itypea(length),itypeb(length)
550 if (itypea(i).ne.itypeb(i)) then
558 c-----------------------------------------------------------------------------
559 subroutine read_bridge
560 C Read information about disulfide bridges.
563 include 'COMMON.IOUNITS'
566 include 'COMMON.INTERACT'
567 include 'COMMON.LOCAL'
568 include 'COMMON.NAMES'
569 include 'COMMON.CHAIN'
570 include 'COMMON.FFIELD'
571 include 'COMMON.SBRIDGE'
572 include 'COMMON.HEADER'
573 include 'COMMON.CONTROL'
574 include 'COMMON.TIME1'
576 include 'COMMON.INFO'
579 C Read bridging residues.
580 read (inp,*) ns,(iss(i),i=1,ns)
582 C Check whether the specified bridging residues are cystines.
584 if (itype(iss(i)).ne.1) then
585 write (iout,'(2a,i3,a)')
586 & 'Do you REALLY think that the residue ',
587 & restyp(itype(iss(i))),i,
588 & ' can form a disulfide bridge?!!!'
589 write (*,'(2a,i3,a)')
590 & 'Do you REALLY think that the residue ',
591 & restyp(itype(iss(i))),i,
592 & ' can form a disulfide bridge?!!!'
594 call mp_stopall(error_msg)
600 C Read preformed bridges.
602 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
605 C Check if the residues involved in bridges are in the specified list of
609 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
610 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
611 write (iout,'(a,i3,a)') 'Disulfide pair',i,
612 & ' contains residues present in other pairs.'
613 write (*,'(a,i3,a)') 'Disulfide pair',i,
614 & ' contains residues present in other pairs.'
616 call mp_stopall(error_msg)
623 if (ihpb(i).eq.iss(j)) goto 10
625 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
628 if (jhpb(i).eq.iss(j)) goto 20
630 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
643 c----------------------------------------------------------------------------
644 subroutine read_angles(kanal,*)
649 include 'COMMON.CHAIN'
650 include 'COMMON.IOUNITS'
652 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
653 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
654 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
655 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
657 theta(i)=deg2rad*theta(i)
658 phi(i)=deg2rad*phi(i)
659 alph(i)=deg2rad*alph(i)
660 omeg(i)=deg2rad*omeg(i)
665 c----------------------------------------------------------------------------
666 subroutine reada(rekord,lancuch,wartosc,default)
668 character*(*) rekord,lancuch
669 double precision wartosc,default
672 iread=index(rekord,lancuch)
677 iread=iread+ilen(lancuch)+1
678 read (rekord(iread:),*) wartosc
681 c----------------------------------------------------------------------------
682 subroutine multreada(rekord,lancuch,tablica,dim,default)
685 double precision tablica(dim),default
686 character*(*) rekord,lancuch
692 iread=index(rekord,lancuch)
693 if (iread.eq.0) return
694 iread=iread+ilen(lancuch)+1
695 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
698 c----------------------------------------------------------------------------
699 subroutine readi(rekord,lancuch,wartosc,default)
701 character*(*) rekord,lancuch
702 integer wartosc,default
705 iread=index(rekord,lancuch)
710 iread=iread+ilen(lancuch)+1
711 read (rekord(iread:),*) wartosc
714 C----------------------------------------------------------------------
715 subroutine multreadi(rekord,lancuch,tablica,dim,default)
718 integer tablica(dim),default
719 character*(*) rekord,lancuch
726 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
727 if (iread.eq.0) return
728 iread=iread+ilen(lancuch)+1
729 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
733 c----------------------------------------------------------------------------
734 subroutine card_concat(card)
736 include 'COMMON.IOUNITS'
738 character*80 karta,ucase
740 read (inp,'(a)') karta
743 do while (karta(80:80).eq.'&')
744 card=card(:ilen(card)+1)//karta(:79)
745 read (inp,'(a)') karta
748 card=card(:ilen(card)+1)//karta
751 c----------------------------------------------------------------------------
760 include 'COMMON.IOUNITS'
761 include 'COMMON.CONTROL'
762 integer lenpre,lenpot,ilen
764 character*16 cformat,cprint
766 integer lenint,lenout
767 call getenv('INPUT',prefix)
768 call getenv('OUTPUT',prefout)
769 call getenv('INTIN',prefintin)
770 call getenv('COORD',cformat)
771 call getenv('PRINTCOOR',cprint)
772 call getenv('SCRATCHDIR',scratchdir)
775 if (index(ucase(cformat),'CX').gt.0) then
782 lenint=ilen(prefintin)
783 C Get the names and open the input files
784 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
786 write (liczba,'(bz,i3.3)') me
787 outname=prefout(:lenout)//'_clust.out_'//liczba
789 outname=prefout(:lenout)//'_clust.out'
792 intinname=prefintin(:lenint)//'.bx'
793 else if (from_cx) then
794 intinname=prefintin(:lenint)//'.cx'
796 intinname=prefintin(:lenint)//'.int'
798 rmsname=prefintin(:lenint)//'.rms'
799 open (jplot,file=prefout(:ilen(prefout))//'.tex',
801 open (jrms,file=rmsname,status='unknown')
802 open(iout,file=outname,status='unknown')
803 C Get parameter filenames and open the parameter files.
804 call getenv('BONDPAR',bondname)
805 open (ibond,file=bondname,status='old')
806 call getenv('THETPAR',thetname)
807 open (ithep,file=thetname,status='old')
808 call getenv('ROTPAR',rotname)
809 open (irotam,file=rotname,status='old')
810 call getenv('TORPAR',torname)
811 open (itorp,file=torname,status='old')
812 call getenv('TORDPAR',tordname)
813 open (itordp,file=tordname,status='old')
814 call getenv('FOURIER',fouriername)
815 open (ifourier,file=fouriername,status='old')
816 call getenv('ELEPAR',elename)
817 open (ielep,file=elename,status='old')
818 call getenv('SIDEPAR',sidename)
819 open (isidep,file=sidename,status='old')
820 call getenv('SIDEP',sidepname)
821 open (isidep1,file=sidepname,status="old")
822 call getenv('SCCORPAR',sccorname)
823 open (isccor,file=sccorname,status="old")
824 call getenv('LIPTRANPAR',liptranname)
825 open (iliptranpar,file=liptranname,status='old')
828 C 8/9/01 In the newest version SCp interaction constants are read from a file
829 C Use -DOLDSCP to use hard-coded constants instead.
831 call getenv('SCPPAR',scpname)
832 open (iscpp,file=scpname,status='old')
836 subroutine read_dist_constr
837 implicit real*8 (a-h,o-z)
842 include 'COMMON.CONTROL'
843 include 'COMMON.CHAIN'
844 include 'COMMON.IOUNITS'
845 include 'COMMON.SBRIDGE'
846 integer ifrag_(2,100),ipair_(2,100)
847 double precision wfrag_(100),wpair_(100)
848 character*500 controlcard
849 logical lprn /.true./
850 write (iout,*) "Calling read_dist_constr"
851 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
853 write(iout,*) "TU sie wywalam?"
854 call card_concat(controlcard)
855 write (iout,*) controlcard
857 call readi(controlcard,"NFRAG",nfrag_,0)
858 call readi(controlcard,"NPAIR",npair_,0)
859 call readi(controlcard,"NDIST",ndist_,0)
860 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
861 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
862 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
863 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
864 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
865 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
866 write (iout,*) "IFRAG"
868 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
870 write (iout,*) "IPAIR"
872 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
876 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
877 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
878 & ifrag_(2,i)=nstart_sup+nsup-1
879 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
881 if (wfrag_(i).gt.0.0d0) then
882 do j=ifrag_(1,i),ifrag_(2,i)-1
884 write (iout,*) "j",j," k",k
886 if (constr_dist.eq.1) then
891 forcon(nhpb)=wfrag_(i)
892 else if (constr_dist.eq.2) then
893 if (ddjk.le.dist_cut) then
898 forcon(nhpb)=wfrag_(i)
905 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
908 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
909 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
915 if (wpair_(i).gt.0.0d0) then
923 do j=ifrag_(1,ii),ifrag_(2,ii)
924 do k=ifrag_(1,jj),ifrag_(2,jj)
928 forcon(nhpb)=wpair_(i)
930 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
931 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
937 if (constr_dist.eq.11) then
938 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
939 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
940 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
941 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
942 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
944 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
946 if (forcon(nhpb+1).gt.0.0d0) then
948 if (ibecarb(i).gt.0) then
952 if (dhpb(nhpb).eq.0.0d0)
953 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
954 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
955 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
956 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
969 c-------------------------------------------------------------------------------
970 subroutine read_saxs_constr
971 implicit real*8 (a-h,o-z)
976 include 'COMMON.CONTROL'
977 include 'COMMON.CHAIN'
978 include 'COMMON.IOUNITS'
979 include 'COMMON.SBRIDGE'
980 double precision cm(3)
982 write (iout,*) "Calling read_saxs nsaxs",nsaxs
984 if (saxs_mode.eq.0) then
985 c SAXS distance distribution
987 read(inp,*) distsaxs(i),Psaxs(i)
991 Cnorm = Cnorm + Psaxs(i)
993 write (iout,*) "Cnorm",Cnorm
995 Psaxs(i)=Psaxs(i)/Cnorm
997 write (iout,*) "Normalized distance distribution from SAXS"
999 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1003 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1005 write (iout,*) "Wsaxs0",Wsaxs0
1009 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1016 cm(j)=cm(j)+Csaxs(j,i)
1024 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1027 write (iout,*) "SAXS sphere coordinates"
1029 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
1034 c====-------------------------------------------------------------------
1035 subroutine read_constr_homology(lprn)
1037 include 'DIMENSIONS'
1041 include 'COMMON.SETUP'
1042 include 'COMMON.CONTROL'
1043 include 'COMMON.CHAIN'
1044 include 'COMMON.IOUNITS'
1045 include 'COMMON.GEO'
1046 include 'COMMON.INTERACT'
1047 include 'COMMON.NAMES'
1048 include 'COMMON.HOMRESTR'
1050 c For new homol impl
1052 include 'COMMON.VAR'
1053 c include 'include_unres/COMMON.VAR'
1056 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1058 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1059 c & sigma_odl_temp(maxres,maxres,max_template)
1061 character*24 model_ki_dist, model_ki_angle
1062 character*500 controlcard
1063 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1064 integer idomain(max_template,maxres)
1068 logical unres_pdb,liiflag
1070 c FP - Nov. 2014 Temporary specifications for new vars
1072 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
1074 double precision, dimension (max_template,maxres) :: rescore
1075 double precision, dimension (max_template,maxres) :: rescore2
1076 double precision, dimension (max_template,maxres) :: rescore3
1077 character*24 tpl_k_rescore
1078 c -----------------------------------------------------------------
1079 c Reading multiple PDB ref structures and calculation of retraints
1080 c not using pre-computed ones stored in files model_ki_{dist,angle}
1082 c -----------------------------------------------------------------
1085 c Alternative: reading from input
1087 write (iout,*) "BEGIN READ HOMOLOGY INFO"
1094 call card_concat(controlcard)
1095 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1096 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1097 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1098 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1099 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1100 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1101 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1102 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1103 if (homol_nset.gt.1)then
1104 call readi(controlcard,"ISET",iset,1)
1105 call card_concat(controlcard)
1106 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1109 waga_homology(1)=1.0
1113 write(iout,*) "read_constr_homology iset",iset
1114 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1121 cd write (iout,*) "nnt",nnt," nct",nct
1131 c Reading HM global scores (prob not required)
1134 do k=1,constr_homology
1138 c open (4,file="HMscore")
1139 c do k=1,constr_homology
1140 c read (4,*,end=521) hmscore_tmp
1141 c hmscore(k)=hmscore_tmp ! Another transformation can be used
1142 c write(*,*) "Model", k, ":", hmscore(k)
1153 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1155 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1156 do k=1,constr_homology
1158 read(inp,'(a)') pdbfile
1159 c write (iout,*) "k ",k," pdbfile ",pdbfile
1160 c Next stament causes error upon compilation (?)
1161 c if(me.eq.king.or. .not. out1file)
1162 c write (iout,'(2a)') 'PDB data will be read from file ',
1163 c & pdbfile(:ilen(pdbfile))
1164 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1165 & pdbfile(:ilen(pdbfile))
1166 open(ipdbin,file=pdbfile,status='old',err=33)
1168 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1169 & pdbfile(:ilen(pdbfile))
1172 c print *,'Begin reading pdb data'
1174 c Files containing res sim or local scores (former containing sigmas)
1177 write(kic2,'(bz,i2.2)') k
1179 tpl_k_rescore="template"//kic2//".sco"
1182 call readpdb_template(k)
1185 crefjlee(j,i)=c(j,i)
1190 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1191 & (crefjlee(j,i+nres),j=1,3)
1193 write (iout,*) "READ HOMOLOGY INFO"
1194 write (iout,*) "read_constr_homology x: after reading pdb file"
1195 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1196 write (iout,*) "waga_dist",waga_dist
1197 write (iout,*) "waga_angle",waga_angle
1198 write (iout,*) "waga_theta",waga_theta
1199 write (iout,*) "waga_d",waga_d
1200 write (iout,*) "dist_cut",dist_cut
1209 c Distance restraints
1212 C Copy the coordinates from reference coordinates (?)
1216 c write (iout,*) "c(",j,i,") =",c(j,i)
1220 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1222 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1223 open (ientin,file=tpl_k_rescore,status='old')
1224 if (nnt.gt.1) rescore(k,1)=0.0d0
1225 do irec=nnt,nct ! loop for reading res sim
1226 if (read2sigma) then
1227 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1228 & rescore3_tmp,idomain_tmp
1230 idomain(k,i_tmp)=idomain_tmp
1231 rescore(k,i_tmp)=rescore_tmp
1232 rescore2(k,i_tmp)=rescore2_tmp
1233 rescore3(k,i_tmp)=rescore3_tmp
1234 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1235 & i_tmp,rescore2_tmp,rescore_tmp,
1236 & rescore3_tmp,idomain_tmp
1239 read (ientin,*,end=1401) rescore_tmp
1241 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1242 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1243 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1248 if (waga_dist.ne.0.0d0) then
1256 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1257 c write (iout,*) k,i,j,distal,dist2_cut
1259 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1260 & .and. distal.le.dist2_cut ) then
1266 c write (iout,*) "k",k
1267 c write (iout,*) "i",i," j",j," constr_homology",
1272 if (read2sigma) then
1275 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1277 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1278 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1279 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1281 if (odl(k,ii).le.dist_cut) then
1282 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1285 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1286 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1288 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1289 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1293 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1296 l_homo(k,ii)=.false.
1303 c Theta, dihedral and SC retraints
1305 if (waga_angle.gt.0.0d0) then
1306 c open (ientin,file=tpl_k_sigma_dih,status='old')
1307 c do irec=1,maxres-3 ! loop for reading sigma_dih
1308 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1309 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1310 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1311 c & sigma_dih(k,i+nnt-1)
1316 if (idomain(k,i).eq.0) then
1320 dih(k,i)=phiref(i) ! right?
1321 c read (ientin,*) sigma_dih(k,i) ! original variant
1322 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1323 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1324 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1325 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1326 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1328 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1329 & rescore(k,i-2)+rescore(k,i-3))/4.0
1330 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1331 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1332 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1333 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1334 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1335 c Instead of res sim other local measure of b/b str reliability possible
1336 if (sigma_dih(k,i).ne.0)
1337 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1338 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1343 if (waga_theta.gt.0.0d0) then
1344 c open (ientin,file=tpl_k_sigma_theta,status='old')
1345 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1346 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1347 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1348 c & sigma_theta(k,i+nnt-1)
1353 do i = nnt+2,nct ! right? without parallel.
1354 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1355 c do i=ithet_start,ithet_end ! with FG parallel.
1356 if (idomain(k,i).eq.0) then
1357 sigma_theta(k,i)=0.0
1360 thetatpl(k,i)=thetaref(i)
1361 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1362 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1363 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1364 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1365 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1366 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1367 & rescore(k,i-2))/3.0
1368 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1369 if (sigma_theta(k,i).ne.0)
1370 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1372 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1373 c rescore(k,i-2) ! right expression ?
1374 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1378 if (waga_d.gt.0.0d0) then
1379 c open (ientin,file=tpl_k_sigma_d,status='old')
1380 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1381 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1382 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1383 c & sigma_d(k,i+nnt-1)
1387 do i = nnt,nct ! right? without parallel.
1388 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1389 c do i=loc_start,loc_end ! with FG parallel.
1390 if (itype(i).eq.10) cycle
1391 if (idomain(k,i).eq.0 ) then
1398 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1399 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1400 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1401 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1402 sigma_d(k,i)=rescore3(k,i) ! right expression ?
1403 if (sigma_d(k,i).ne.0)
1404 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1406 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1407 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1408 c read (ientin,*) sigma_d(k,i) ! 1st variant
1413 c remove distance restraints not used in any model from the list
1414 c shift data in all arrays
1416 if (waga_dist.ne.0.0d0) then
1422 if (ii_in_use(ii).eq.0.and.liiflag) then
1426 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1427 & .not.liiflag.and.ii.eq.lim_odl) then
1428 if (ii.eq.lim_odl) then
1429 iishift=ii-iistart+1
1434 do ki=iistart,lim_odl-iishift
1435 ires_homo(ki)=ires_homo(ki+iishift)
1436 jres_homo(ki)=jres_homo(ki+iishift)
1437 ii_in_use(ki)=ii_in_use(ki+iishift)
1438 do k=1,constr_homology
1439 odl(k,ki)=odl(k,ki+iishift)
1440 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1441 l_homo(k,ki)=l_homo(k,ki+iishift)
1445 lim_odl=lim_odl-iishift
1450 if (constr_homology.gt.0) call homology_partition
1451 if (constr_homology.gt.0) call init_int_table
1452 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1453 cd & "lim_xx=",lim_xx
1454 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1455 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1459 if (.not.lprn) return
1460 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1461 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1462 write (iout,*) "Distance restraints from templates"
1464 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1465 & ii,ires_homo(ii),jres_homo(ii),
1466 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1467 & ki=1,constr_homology)
1469 write (iout,*) "Dihedral angle restraints from templates"
1471 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1472 & (rad2deg*dih(ki,i),
1473 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1475 write (iout,*) "Virtual-bond angle restraints from templates"
1477 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1478 & (rad2deg*thetatpl(ki,i),
1479 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1481 write (iout,*) "SC restraints from templates"
1483 write(iout,'(i5,100(4f8.2,4x))') i,
1484 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1485 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1488 c -----------------------------------------------------------------