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 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
138 c--------------------------------------------------------------------------
141 C Read molecular data.
145 include 'COMMON.IOUNITS'
148 include 'COMMON.INTERACT'
149 include 'COMMON.LOCAL'
150 include 'COMMON.NAMES'
151 include 'COMMON.CHAIN'
152 include 'COMMON.FFIELD'
153 include 'COMMON.SBRIDGE'
154 include 'COMMON.HEADER'
155 include 'COMMON.CONTROL'
156 include 'COMMON.CONTACTS'
157 include 'COMMON.TIME1'
158 include 'COMMON.TORCNSTR'
159 include 'COMMON.SHIELD'
161 include 'COMMON.INFO'
163 character*4 sequence(maxres)
164 character*800 weightcard
166 double precision x(maxvar)
167 integer itype_pdb(maxres)
169 integer i,j,kkk,i1,i2,it1,it2
173 C Read weights of the subsequent energy terms.
174 call card_concat(weightcard)
175 write(iout,*) weightcard
176 call reada(weightcard,'WSC',wsc,1.0d0)
177 call reada(weightcard,'WLONG',wsc,wsc)
178 call reada(weightcard,'WSCP',wscp,1.0d0)
179 call reada(weightcard,'WELEC',welec,1.0D0)
180 call reada(weightcard,'WVDWPP',wvdwpp,welec)
181 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
182 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
183 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
184 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
185 call reada(weightcard,'WTURN3',wturn3,1.0D0)
186 call reada(weightcard,'WTURN4',wturn4,1.0D0)
187 call reada(weightcard,'WTURN6',wturn6,1.0D0)
188 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
189 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
190 call reada(weightcard,'WBOND',wbond,1.0D0)
191 call reada(weightcard,'WTOR',wtor,1.0D0)
192 call reada(weightcard,'WTORD',wtor_d,1.0D0)
193 call reada(weightcard,'WANG',wang,1.0D0)
194 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
195 call reada(weightcard,'WSAXS',wsaxs,0.0D0)
196 call reada(weightcard,'SCAL14',scal14,0.4D0)
197 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
198 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
199 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
200 if (index(weightcard,'SOFT').gt.0) ipot=6
201 call reada(weightcard,"D0CM",d0cm,3.78d0)
202 call reada(weightcard,"AKCM",akcm,15.1d0)
203 call reada(weightcard,"AKTH",akth,11.0d0)
204 call reada(weightcard,"AKCT",akct,12.0d0)
205 call reada(weightcard,"V1SS",v1ss,-1.08d0)
206 call reada(weightcard,"V2SS",v2ss,7.61d0)
207 call reada(weightcard,"V3SS",v3ss,13.7d0)
208 call reada(weightcard,"EBR",ebr,-5.50D0)
209 call reada(weightcard,'WSHIELD',wshield,1.0d0)
210 write(iout,*) 'WSHIELD',wshield
211 call reada(weightcard,'WLT',wliptran,0.0D0)
212 call reada(weightcard,"ATRISS",atriss,0.301D0)
213 call reada(weightcard,"BTRISS",btriss,0.021D0)
214 call reada(weightcard,"CTRISS",ctriss,1.001D0)
215 call reada(weightcard,"DTRISS",dtriss,1.001D0)
216 write (iout,*) "ATRISS=", atriss
217 write (iout,*) "BTRISS=", btriss
218 write (iout,*) "CTRISS=", ctriss
219 write (iout,*) "DTRISS=", dtriss
220 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
222 dyn_ss_mask(i)=.false.
226 dyn_ssbond_ij(i,j)=1.0d300
229 call reada(weightcard,"HT",Ht,0.0D0)
231 ss_depth=ebr/wsc-0.25*eps(1,1)
232 Ht=Ht/wsc-0.25*eps(1,1)
233 akcm=akcm*wstrain/wsc
234 akth=akth*wstrain/wsc
235 akct=akct*wstrain/wsc
236 v1ss=v1ss*wstrain/wsc
237 v2ss=v2ss*wstrain/wsc
238 v3ss=v3ss*wstrain/wsc
240 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
242 write (iout,'(/a)') "Disulfide bridge parameters:"
243 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
244 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
245 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
246 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
249 C 12/1/95 Added weight for the multi-body term WCORR
250 call reada(weightcard,'WCORRH',wcorr,1.0D0)
251 if (wcorr4.gt.0.0d0) wcorr=wcorr4
271 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
272 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
273 & wturn4,wturn6,wsccor
274 10 format (/'Energy-term weights (unscaled):'//
275 & 'WSCC= ',f10.6,' (SC-SC)'/
276 & 'WSCP= ',f10.6,' (SC-p)'/
277 & 'WELEC= ',f10.6,' (p-p electr)'/
278 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
279 & 'WBOND= ',f10.6,' (stretching)'/
280 & 'WANG= ',f10.6,' (bending)'/
281 & 'WSCLOC= ',f10.6,' (SC local)'/
282 & 'WTOR= ',f10.6,' (torsional)'/
283 & 'WTORD= ',f10.6,' (double torsional)'/
284 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
285 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
286 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
287 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
288 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
289 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
290 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
291 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
292 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
294 if (wcorr4.gt.0.0d0) then
295 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
296 & 'between contact pairs of peptide groups'
297 write (iout,'(2(a,f5.3/))')
298 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
299 & 'Range of quenching the correlation terms:',2*delt_corr
300 else if (wcorr.gt.0.0d0) then
301 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
302 & 'between contact pairs of peptide groups'
304 write (iout,'(a,f8.3)')
305 & 'Scaling factor of 1,4 SC-p interactions:',scal14
306 write (iout,'(a,f8.3)')
307 & 'General scaling factor of SC-p interactions:',scalscp
308 r0_corr=cutoff_corr-delt_corr
310 aad(i,1)=scalscp*aad(i,1)
311 aad(i,2)=scalscp*aad(i,2)
312 bad(i,1)=scalscp*bad(i,1)
313 bad(i,2)=scalscp*bad(i,2)
317 print *,'indpdb=',indpdb,' pdbref=',pdbref
319 C Read sequence if not taken from the pdb file.
320 if (iscode.gt.0) then
321 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
323 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
325 C Convert sequence to numeric code
327 itype(i)=rescode(i,sequence(i),iscode)
330 print '(20i4)',(itype(i),i=1,nres)
334 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
336 if (itype(i).eq.ntyp1) then
340 else if (iabs(itype(i+1)).ne.20) then
342 else if (iabs(itype(i)).ne.20) then
349 write (iout,*) "ITEL"
351 write (iout,*) i,itype(i),itel(i)
354 print *,'Call Read_Bridge.'
356 C this fragment reads diheadral constrains
357 if (with_dihed_constr) then
359 read (inp,*) ndih_constr
360 if (ndih_constr.gt.0) then
362 C write (iout,*) 'FTORS',ftors
363 C ftors is the force constant for torsional quartic constrains
364 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
367 & 'There are',ndih_constr,' constraints on phi angles.'
369 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
373 phi0(i)=deg2rad*phi0(i)
374 drange(i)=deg2rad*drange(i)
376 endif ! endif ndif_constr.gt.0
377 endif ! with_dihed_constr
378 if (with_theta_constr) then
379 C with_theta_constr is keyword allowing for occurance of theta constrains
380 read (inp,*) ntheta_constr
381 C ntheta_constr is the number of theta constrains
382 if (ntheta_constr.gt.0) then
384 read (inp,*) (itheta_constr(i),theta_constr0(i),
385 & theta_drange(i),for_thet_constr(i),
387 C the above code reads from 1 to ntheta_constr
388 C itheta_constr(i) residue i for which is theta_constr
389 C theta_constr0 the global minimum value
390 C theta_drange is range for which there is no energy penalty
391 C for_thet_constr is the force constant for quartic energy penalty
393 C if(me.eq.king.or..not.out1file)then
395 & 'There are',ntheta_constr,' constraints on phi angles.'
397 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
403 theta_constr0(i)=deg2rad*theta_constr0(i)
404 theta_drange(i)=deg2rad*theta_drange(i)
406 C if(me.eq.king.or..not.out1file)
407 C & write (iout,*) 'FTORS',ftors
408 C do i=1,ntheta_constr
409 C ii = itheta_constr(i)
410 C thetabound(1,ii) = phi0(i)-drange(i)
411 C thetabound(2,ii) = phi0(i)+drange(i)
413 endif ! ntheta_constr.gt.0
414 endif! with_theta_constr
418 print *,'NNT=',NNT,' NCT=',NCT
419 if (itype(1).eq.ntyp1) nnt=2
420 if (itype(nres).eq.ntyp1) nct=nct-1
421 if (nstart.lt.nnt) nstart=nnt
422 if (nend.gt.nct .or. nend.eq.0) nend=nct
423 write (iout,*) "nstart",nstart," nend",nend
424 write (iout,*) "calling read_saxs_consrtr",nsaxs
425 if (nsaxs.gt.0) call read_saxs_constr
427 if (constr_homology.gt.0) then
428 call read_constr_homology(print_homology_restraints)
432 c read(inp,'(a)') pdbfile
433 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
434 c open(ipdbin,file=pdbfile,status='old',err=33)
436 c 33 write (iout,'(a)') 'Error opening PDB file.'
439 c print *,'Begin reading pdb data'
441 c print *,'Finished reading pdb data'
442 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
444 c itype_pdb(i)=itype(i)
447 c write (iout,'(a,i3)') 'nsup=',nsup
449 c if (nsup.le.(nct-nnt+1)) then
450 c do i=0,nct-nnt+1-nsup
451 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
457 c & 'Error - sequences to be superposed do not match.'
460 c do i=0,nsup-(nct-nnt+1)
461 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
463 c nstart_sup=nstart_sup+i
469 c & 'Error - sequences to be superposed do not match.'
472 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
473 c & ' nstart_seq=',nstart_seq
477 write (iout,*) "molread: REFSTR",refstr
479 if (.not.pdbref) then
480 call read_angles(inp,*38)
482 38 write (iout,'(a)') 'Error reading reference structure.'
484 call mp_stopall(Error_Msg)
486 stop 'Error reading reference structure'
499 call contact(.true.,ncont_ref,icont_ref)
502 C write (iout,'(/a,i3,a)')
503 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
504 write (iout,'(20i4)') (iss(i),i=1,ns)
506 write(iout,*)"Running with dynamic disulfide-bond formation"
508 write (iout,'(/a/)') 'Pre-formed links are:'
514 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
515 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
521 if (ns.gt.0.and.dyn_ss) then
525 forcon(i-nss)=forcon(i)
532 dyn_ss_mask(iss(i))=.true.
535 c Read distance restraints
536 if (constr_dist.gt.0) then
537 call read_dist_constr
542 c-----------------------------------------------------------------------------
543 logical function seq_comp(itypea,itypeb,length)
545 integer length,itypea(length),itypeb(length)
548 if (itypea(i).ne.itypeb(i)) then
556 c-----------------------------------------------------------------------------
557 subroutine read_bridge
558 C Read information about disulfide bridges.
561 include 'COMMON.IOUNITS'
564 include 'COMMON.INTERACT'
565 include 'COMMON.LOCAL'
566 include 'COMMON.NAMES'
567 include 'COMMON.CHAIN'
568 include 'COMMON.FFIELD'
569 include 'COMMON.SBRIDGE'
570 include 'COMMON.HEADER'
571 include 'COMMON.CONTROL'
572 include 'COMMON.TIME1'
574 include 'COMMON.INFO'
577 C Read bridging residues.
578 read (inp,*) ns,(iss(i),i=1,ns)
580 C Check whether the specified bridging residues are cystines.
582 if (itype(iss(i)).ne.1) then
583 write (iout,'(2a,i3,a)')
584 & 'Do you REALLY think that the residue ',
585 & restyp(itype(iss(i))),i,
586 & ' can form a disulfide bridge?!!!'
587 write (*,'(2a,i3,a)')
588 & 'Do you REALLY think that the residue ',
589 & restyp(itype(iss(i))),i,
590 & ' can form a disulfide bridge?!!!'
592 call mp_stopall(error_msg)
598 C Read preformed bridges.
600 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
603 C Check if the residues involved in bridges are in the specified list of
607 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
608 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
609 write (iout,'(a,i3,a)') 'Disulfide pair',i,
610 & ' contains residues present in other pairs.'
611 write (*,'(a,i3,a)') 'Disulfide pair',i,
612 & ' contains residues present in other pairs.'
614 call mp_stopall(error_msg)
621 if (ihpb(i).eq.iss(j)) goto 10
623 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
626 if (jhpb(i).eq.iss(j)) goto 20
628 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
641 c----------------------------------------------------------------------------
642 subroutine read_angles(kanal,*)
647 include 'COMMON.CHAIN'
648 include 'COMMON.IOUNITS'
650 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
651 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
652 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
653 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
655 theta(i)=deg2rad*theta(i)
656 phi(i)=deg2rad*phi(i)
657 alph(i)=deg2rad*alph(i)
658 omeg(i)=deg2rad*omeg(i)
663 c----------------------------------------------------------------------------
664 subroutine reada(rekord,lancuch,wartosc,default)
666 character*(*) rekord,lancuch
667 double precision wartosc,default
670 iread=index(rekord,lancuch)
675 iread=iread+ilen(lancuch)+1
676 read (rekord(iread:),*) wartosc
679 c----------------------------------------------------------------------------
680 subroutine multreada(rekord,lancuch,tablica,dim,default)
683 double precision tablica(dim),default
684 character*(*) rekord,lancuch
690 iread=index(rekord,lancuch)
691 if (iread.eq.0) return
692 iread=iread+ilen(lancuch)+1
693 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
696 c----------------------------------------------------------------------------
697 subroutine readi(rekord,lancuch,wartosc,default)
699 character*(*) rekord,lancuch
700 integer wartosc,default
703 iread=index(rekord,lancuch)
708 iread=iread+ilen(lancuch)+1
709 read (rekord(iread:),*) wartosc
712 C----------------------------------------------------------------------
713 subroutine multreadi(rekord,lancuch,tablica,dim,default)
716 integer tablica(dim),default
717 character*(*) rekord,lancuch
724 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
725 if (iread.eq.0) return
726 iread=iread+ilen(lancuch)+1
727 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
731 c----------------------------------------------------------------------------
732 subroutine card_concat(card)
734 include 'COMMON.IOUNITS'
736 character*80 karta,ucase
738 read (inp,'(a)') karta
741 do while (karta(80:80).eq.'&')
742 card=card(:ilen(card)+1)//karta(:79)
743 read (inp,'(a)') karta
746 card=card(:ilen(card)+1)//karta
749 c----------------------------------------------------------------------------
758 include 'COMMON.IOUNITS'
759 include 'COMMON.CONTROL'
760 integer lenpre,lenpot,ilen
762 character*16 cformat,cprint
764 integer lenint,lenout
765 call getenv('INPUT',prefix)
766 call getenv('OUTPUT',prefout)
767 call getenv('INTIN',prefintin)
768 call getenv('COORD',cformat)
769 call getenv('PRINTCOOR',cprint)
770 call getenv('SCRATCHDIR',scratchdir)
773 if (index(ucase(cformat),'CX').gt.0) then
780 lenint=ilen(prefintin)
781 C Get the names and open the input files
782 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
784 write (liczba,'(bz,i3.3)') me
785 outname=prefout(:lenout)//'_clust.out_'//liczba
787 outname=prefout(:lenout)//'_clust.out'
790 intinname=prefintin(:lenint)//'.bx'
791 else if (from_cx) then
792 intinname=prefintin(:lenint)//'.cx'
794 intinname=prefintin(:lenint)//'.int'
796 rmsname=prefintin(:lenint)//'.rms'
797 open (jplot,file=prefout(:ilen(prefout))//'.tex',
799 open (jrms,file=rmsname,status='unknown')
800 open(iout,file=outname,status='unknown')
801 C Get parameter filenames and open the parameter files.
802 call getenv('BONDPAR',bondname)
803 open (ibond,file=bondname,status='old')
804 call getenv('THETPAR',thetname)
805 open (ithep,file=thetname,status='old')
806 call getenv('ROTPAR',rotname)
807 open (irotam,file=rotname,status='old')
808 call getenv('TORPAR',torname)
809 open (itorp,file=torname,status='old')
810 call getenv('TORDPAR',tordname)
811 open (itordp,file=tordname,status='old')
812 call getenv('FOURIER',fouriername)
813 open (ifourier,file=fouriername,status='old')
814 call getenv('ELEPAR',elename)
815 open (ielep,file=elename,status='old')
816 call getenv('SIDEPAR',sidename)
817 open (isidep,file=sidename,status='old')
818 call getenv('SIDEP',sidepname)
819 open (isidep1,file=sidepname,status="old")
820 call getenv('SCCORPAR',sccorname)
821 open (isccor,file=sccorname,status="old")
822 call getenv('LIPTRANPAR',liptranname)
823 open (iliptranpar,file=liptranname,status='old')
826 C 8/9/01 In the newest version SCp interaction constants are read from a file
827 C Use -DOLDSCP to use hard-coded constants instead.
829 call getenv('SCPPAR',scpname)
830 open (iscpp,file=scpname,status='old')
834 subroutine read_dist_constr
835 implicit real*8 (a-h,o-z)
840 include 'COMMON.CONTROL'
841 include 'COMMON.CHAIN'
842 include 'COMMON.IOUNITS'
843 include 'COMMON.SBRIDGE'
844 integer ifrag_(2,100),ipair_(2,100)
845 double precision wfrag_(100),wpair_(100)
846 character*500 controlcard
847 logical lprn /.true./
848 write (iout,*) "Calling read_dist_constr"
849 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
851 write(iout,*) "TU sie wywalam?"
852 call card_concat(controlcard)
853 write (iout,*) controlcard
855 call readi(controlcard,"NFRAG",nfrag_,0)
856 call readi(controlcard,"NPAIR",npair_,0)
857 call readi(controlcard,"NDIST",ndist_,0)
858 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
859 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
860 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
861 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
862 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
863 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
864 write (iout,*) "IFRAG"
866 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
868 write (iout,*) "IPAIR"
870 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
874 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
875 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
876 & ifrag_(2,i)=nstart_sup+nsup-1
877 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
879 if (wfrag_(i).gt.0.0d0) then
880 do j=ifrag_(1,i),ifrag_(2,i)-1
882 write (iout,*) "j",j," k",k
884 if (constr_dist.eq.1) then
889 forcon(nhpb)=wfrag_(i)
890 else if (constr_dist.eq.2) then
891 if (ddjk.le.dist_cut) then
896 forcon(nhpb)=wfrag_(i)
903 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
906 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
907 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
913 if (wpair_(i).gt.0.0d0) then
921 do j=ifrag_(1,ii),ifrag_(2,ii)
922 do k=ifrag_(1,jj),ifrag_(2,jj)
926 forcon(nhpb)=wpair_(i)
928 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
929 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
935 if (constr_dist.eq.11) then
936 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
937 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
938 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
939 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
940 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
942 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
944 if (forcon(nhpb+1).gt.0.0d0) then
946 if (ibecarb(i).gt.0) then
950 if (dhpb(nhpb).eq.0.0d0)
951 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
952 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
953 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
954 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
967 c-------------------------------------------------------------------------------
968 subroutine read_saxs_constr
969 implicit real*8 (a-h,o-z)
974 include 'COMMON.CONTROL'
975 include 'COMMON.CHAIN'
976 include 'COMMON.IOUNITS'
977 include 'COMMON.SBRIDGE'
978 double precision cm(3)
980 write (iout,*) "Calling read_saxs nsaxs",nsaxs
982 if (saxs_mode.eq.0) then
983 c SAXS distance distribution
985 read(inp,*) distsaxs(i),Psaxs(i)
989 Cnorm = Cnorm + Psaxs(i)
991 write (iout,*) "Cnorm",Cnorm
993 Psaxs(i)=Psaxs(i)/Cnorm
995 write (iout,*) "Normalized distance distribution from SAXS"
997 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1002 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1009 cm(j)=cm(j)+Csaxs(j,i)
1017 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1020 write (iout,*) "SAXS sphere coordinates"
1022 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
1027 c====-------------------------------------------------------------------
1028 subroutine read_constr_homology(lprn)
1030 include 'DIMENSIONS'
1034 include 'COMMON.SETUP'
1035 include 'COMMON.CONTROL'
1036 include 'COMMON.CHAIN'
1037 include 'COMMON.IOUNITS'
1038 include 'COMMON.GEO'
1039 include 'COMMON.INTERACT'
1040 include 'COMMON.NAMES'
1041 include 'COMMON.HOMRESTR'
1043 c For new homol impl
1045 include 'COMMON.VAR'
1046 c include 'include_unres/COMMON.VAR'
1049 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1051 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1052 c & sigma_odl_temp(maxres,maxres,max_template)
1054 character*24 model_ki_dist, model_ki_angle
1055 character*500 controlcard
1056 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1057 integer idomain(max_template,maxres)
1061 logical unres_pdb,liiflag
1063 c FP - Nov. 2014 Temporary specifications for new vars
1065 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
1067 double precision, dimension (max_template,maxres) :: rescore
1068 double precision, dimension (max_template,maxres) :: rescore2
1069 double precision, dimension (max_template,maxres) :: rescore3
1070 character*24 tpl_k_rescore
1071 c -----------------------------------------------------------------
1072 c Reading multiple PDB ref structures and calculation of retraints
1073 c not using pre-computed ones stored in files model_ki_{dist,angle}
1075 c -----------------------------------------------------------------
1078 c Alternative: reading from input
1080 write (iout,*) "BEGIN READ HOMOLOGY INFO"
1087 call card_concat(controlcard)
1088 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1089 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1090 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1091 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1092 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1093 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1094 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1095 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1096 if (homol_nset.gt.1)then
1097 call readi(controlcard,"ISET",iset,1)
1098 call card_concat(controlcard)
1099 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1102 waga_homology(1)=1.0
1106 write(iout,*) "read_constr_homology iset",iset
1107 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1114 cd write (iout,*) "nnt",nnt," nct",nct
1124 c Reading HM global scores (prob not required)
1127 do k=1,constr_homology
1131 c open (4,file="HMscore")
1132 c do k=1,constr_homology
1133 c read (4,*,end=521) hmscore_tmp
1134 c hmscore(k)=hmscore_tmp ! Another transformation can be used
1135 c write(*,*) "Model", k, ":", hmscore(k)
1146 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1148 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1149 do k=1,constr_homology
1151 read(inp,'(a)') pdbfile
1152 c write (iout,*) "k ",k," pdbfile ",pdbfile
1153 c Next stament causes error upon compilation (?)
1154 c if(me.eq.king.or. .not. out1file)
1155 c write (iout,'(2a)') 'PDB data will be read from file ',
1156 c & pdbfile(:ilen(pdbfile))
1157 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1158 & pdbfile(:ilen(pdbfile))
1159 open(ipdbin,file=pdbfile,status='old',err=33)
1161 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1162 & pdbfile(:ilen(pdbfile))
1165 c print *,'Begin reading pdb data'
1167 c Files containing res sim or local scores (former containing sigmas)
1170 write(kic2,'(bz,i2.2)') k
1172 tpl_k_rescore="template"//kic2//".sco"
1175 call readpdb_template(k)
1178 crefjlee(j,i)=c(j,i)
1183 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1184 & (crefjlee(j,i+nres),j=1,3)
1186 write (iout,*) "READ HOMOLOGY INFO"
1187 write (iout,*) "read_constr_homology x: after reading pdb file"
1188 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1189 write (iout,*) "waga_dist",waga_dist
1190 write (iout,*) "waga_angle",waga_angle
1191 write (iout,*) "waga_theta",waga_theta
1192 write (iout,*) "waga_d",waga_d
1193 write (iout,*) "dist_cut",dist_cut
1202 c Distance restraints
1205 C Copy the coordinates from reference coordinates (?)
1209 c write (iout,*) "c(",j,i,") =",c(j,i)
1213 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1215 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1216 open (ientin,file=tpl_k_rescore,status='old')
1217 if (nnt.gt.1) rescore(k,1)=0.0d0
1218 do irec=nnt,nct ! loop for reading res sim
1219 if (read2sigma) then
1220 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1221 & rescore3_tmp,idomain_tmp
1223 idomain(k,i_tmp)=idomain_tmp
1224 rescore(k,i_tmp)=rescore_tmp
1225 rescore2(k,i_tmp)=rescore2_tmp
1226 rescore3(k,i_tmp)=rescore3_tmp
1227 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1228 & i_tmp,rescore2_tmp,rescore_tmp,
1229 & rescore3_tmp,idomain_tmp
1232 read (ientin,*,end=1401) rescore_tmp
1234 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1235 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1236 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1241 if (waga_dist.ne.0.0d0) then
1249 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1250 c write (iout,*) k,i,j,distal,dist2_cut
1252 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1253 & .and. distal.le.dist2_cut ) then
1259 c write (iout,*) "k",k
1260 c write (iout,*) "i",i," j",j," constr_homology",
1265 if (read2sigma) then
1268 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1270 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1271 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1272 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1274 if (odl(k,ii).le.dist_cut) then
1275 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1278 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1279 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1281 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1282 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1286 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1289 l_homo(k,ii)=.false.
1296 c Theta, dihedral and SC retraints
1298 if (waga_angle.gt.0.0d0) then
1299 c open (ientin,file=tpl_k_sigma_dih,status='old')
1300 c do irec=1,maxres-3 ! loop for reading sigma_dih
1301 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1302 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1303 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1304 c & sigma_dih(k,i+nnt-1)
1309 if (idomain(k,i).eq.0) then
1313 dih(k,i)=phiref(i) ! right?
1314 c read (ientin,*) sigma_dih(k,i) ! original variant
1315 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1316 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1317 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1318 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1319 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1321 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1322 & rescore(k,i-2)+rescore(k,i-3))/4.0
1323 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1324 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1325 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1326 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1327 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1328 c Instead of res sim other local measure of b/b str reliability possible
1329 if (sigma_dih(k,i).ne.0)
1330 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1331 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1336 if (waga_theta.gt.0.0d0) then
1337 c open (ientin,file=tpl_k_sigma_theta,status='old')
1338 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1339 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1340 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1341 c & sigma_theta(k,i+nnt-1)
1346 do i = nnt+2,nct ! right? without parallel.
1347 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1348 c do i=ithet_start,ithet_end ! with FG parallel.
1349 if (idomain(k,i).eq.0) then
1350 sigma_theta(k,i)=0.0
1353 thetatpl(k,i)=thetaref(i)
1354 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1355 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1356 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1357 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1358 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1359 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1360 & rescore(k,i-2))/3.0
1361 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1362 if (sigma_theta(k,i).ne.0)
1363 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1365 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1366 c rescore(k,i-2) ! right expression ?
1367 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1371 if (waga_d.gt.0.0d0) then
1372 c open (ientin,file=tpl_k_sigma_d,status='old')
1373 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1374 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1375 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1376 c & sigma_d(k,i+nnt-1)
1380 do i = nnt,nct ! right? without parallel.
1381 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1382 c do i=loc_start,loc_end ! with FG parallel.
1383 if (itype(i).eq.10) cycle
1384 if (idomain(k,i).eq.0 ) then
1391 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1392 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1393 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1394 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1395 sigma_d(k,i)=rescore3(k,i) ! right expression ?
1396 if (sigma_d(k,i).ne.0)
1397 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1399 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1400 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1401 c read (ientin,*) sigma_d(k,i) ! 1st variant
1406 c remove distance restraints not used in any model from the list
1407 c shift data in all arrays
1409 if (waga_dist.ne.0.0d0) then
1415 if (ii_in_use(ii).eq.0.and.liiflag) then
1419 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1420 & .not.liiflag.and.ii.eq.lim_odl) then
1421 if (ii.eq.lim_odl) then
1422 iishift=ii-iistart+1
1427 do ki=iistart,lim_odl-iishift
1428 ires_homo(ki)=ires_homo(ki+iishift)
1429 jres_homo(ki)=jres_homo(ki+iishift)
1430 ii_in_use(ki)=ii_in_use(ki+iishift)
1431 do k=1,constr_homology
1432 odl(k,ki)=odl(k,ki+iishift)
1433 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1434 l_homo(k,ki)=l_homo(k,ki+iishift)
1438 lim_odl=lim_odl-iishift
1443 if (constr_homology.gt.0) call homology_partition
1444 if (constr_homology.gt.0) call init_int_table
1445 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1446 cd & "lim_xx=",lim_xx
1447 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1448 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1452 if (.not.lprn) return
1453 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1454 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1455 write (iout,*) "Distance restraints from templates"
1457 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1458 & ii,ires_homo(ii),jres_homo(ii),
1459 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1460 & ki=1,constr_homology)
1462 write (iout,*) "Dihedral angle restraints from templates"
1464 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1465 & (rad2deg*dih(ki,i),
1466 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1468 write (iout,*) "Virtual-bond angle restraints from templates"
1470 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1471 & (rad2deg*thetatpl(ki,i),
1472 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1474 write (iout,*) "SC restraints from templates"
1476 write(iout,'(i5,100(4f8.2,4x))') i,
1477 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1478 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1481 c -----------------------------------------------------------------