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'
17 character*320 controlcard,ucase
23 read (INP,'(a80)') titel
24 call card_concat(controlcard)
26 call readi(controlcard,'NRES',nres,0)
27 write (iout,*) "NRES",NRES
28 call readi(controlcard,'RESCALE',rescale_mode,2)
29 call readi(controlcard,'PDBOUT',outpdb,0)
30 call readi(controlcard,'MOL2OUT',outmol2,0)
31 refstr=(index(controlcard,'REFSTR').gt.0)
32 write (iout,*) "REFSTR",refstr
33 pdbref=(index(controlcard,'PDBREF').gt.0)
34 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
35 iscode=index(controlcard,'ONE_LETTER')
36 tree=(index(controlcard,'MAKE_TREE').gt.0)
37 min_var=(index(controlcard,'MINVAR').gt.0)
38 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
39 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
40 call readi(controlcard,'NCUT',ncut,0)
41 call readi(controlcard,'NCLUST',nclust,5)
42 call readi(controlcard,'NSTART',nstart,0)
43 call readi(controlcard,'NEND',nend,0)
44 call reada(controlcard,'ECUT',ecut,10.0d0)
45 call reada(controlcard,'PROB',prob_limit,0.99d0)
46 write (iout,*) "Probability limit",prob_limit
47 lgrp=(index(controlcard,'LGRP').gt.0)
48 caonly=(index(controlcard,'CA_ONLY').gt.0)
49 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
51 & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
52 call readi(controlcard,'IOPT',iopt,2)
53 lside = index(controlcard,"SIDE").gt.0
54 efree = index(controlcard,"EFREE").gt.0
55 call readi(controlcard,'NTEMP',nT,1)
56 write (iout,*) "nT",nT
57 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
58 write (iout,*) "nT",nT
59 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
61 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
63 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
64 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
65 lprint_int=index(controlcard,"PRINT_INT") .gt.0
66 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
67 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
68 write (iout,*) "with_dihed_constr ",with_dihed_constr,
69 & " CONSTR_DIST",constr_dist
71 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
72 write (iout,*) "with_homology_constr ",with_dihed_constr,
73 & " CONSTR_HOMOLOGY",constr_homology
74 print_homology_restraints=
75 & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
76 print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
77 print_homology_models=
78 & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
79 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
89 c--------------------------------------------------------------------------
92 C Read molecular data.
96 include 'COMMON.IOUNITS'
99 include 'COMMON.INTERACT'
100 include 'COMMON.LOCAL'
101 include 'COMMON.NAMES'
102 include 'COMMON.CHAIN'
103 include 'COMMON.FFIELD'
104 include 'COMMON.SBRIDGE'
105 include 'COMMON.HEADER'
106 include 'COMMON.CONTROL'
107 include 'COMMON.CONTACTS'
108 include 'COMMON.TIME1'
109 include 'COMMON.TORCNSTR'
111 include 'COMMON.INFO'
113 character*4 sequence(maxres)
114 character*800 weightcard
116 double precision x(maxvar)
117 integer itype_pdb(maxres)
120 write (iout,*) " MOLREAD: NRES",NRES
124 C Read weights of the subsequent energy terms.
125 call card_concat(weightcard)
126 call reada(weightcard,'WSC',wsc,1.0d0)
127 call reada(weightcard,'WLONG',wsc,wsc)
128 call reada(weightcard,'WSCP',wscp,1.0d0)
129 call reada(weightcard,'WELEC',welec,1.0D0)
130 call reada(weightcard,'WVDWPP',wvdwpp,welec)
131 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
132 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
133 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
134 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
135 call reada(weightcard,'WTURN3',wturn3,1.0D0)
136 call reada(weightcard,'WTURN4',wturn4,1.0D0)
137 call reada(weightcard,'WTURN6',wturn6,1.0D0)
138 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
139 call reada(weightcard,'WBOND',wbond,1.0D0)
140 call reada(weightcard,'WTOR',wtor,1.0D0)
141 call reada(weightcard,'WTORD',wtor_d,1.0D0)
142 call reada(weightcard,'WANG',wang,1.0D0)
143 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
144 call reada(weightcard,'SCAL14',scal14,0.4D0)
145 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
146 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
147 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
148 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
149 call reada(weightcard,"D0CM",d0cm,3.78d0)
150 call reada(weightcard,"AKCM",akcm,15.1d0)
151 call reada(weightcard,"AKTH",akth,11.0d0)
152 call reada(weightcard,"AKCT",akct,12.0d0)
153 call reada(weightcard,"V1SS",v1ss,-1.08d0)
154 call reada(weightcard,"V2SS",v2ss,7.61d0)
155 call reada(weightcard,"V3SS",v3ss,13.7d0)
156 call reada(weightcard,"EBR",ebr,-5.50D0)
158 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
159 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
160 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
161 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
162 if (index(weightcard,'SOFT').gt.0) ipot=6
163 C 12/1/95 Added weight for the multi-body term WCORR
164 call reada(weightcard,'WCORRH',wcorr,1.0D0)
167 dyn_ssbond_ij(i,j)=1.0d300
170 call reada(weightcard,"HT",Ht,0.0D0)
172 ss_depth=ebr/wsc-0.25*eps(1,1)
173 Ht=Ht/wsc-0.25*eps(1,1)
174 akcm=akcm*wstrain/wsc
175 akth=akth*wstrain/wsc
176 akct=akct*wstrain/wsc
177 v1ss=v1ss*wstrain/wsc
178 v2ss=v2ss*wstrain/wsc
179 v3ss=v3ss*wstrain/wsc
181 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
183 write (iout,'(/a)') "Disulfide bridge parameters:"
184 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
185 write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
186 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
187 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
188 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
190 write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
191 if (wcorr4.gt.0.0d0) wcorr=wcorr4
210 weights(22)=wdfa_dist
213 weights(25)=wdfa_beta
214 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
215 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
216 & wturn4,wturn6,wsccor,wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
217 10 format (/'Energy-term weights (unscaled):'//
218 & 'WSCC= ',f10.6,' (SC-SC)'/
219 & 'WSCP= ',f10.6,' (SC-p)'/
220 & 'WELEC= ',f10.6,' (p-p electr)'/
221 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
222 & 'WBOND= ',f10.6,' (stretching)'/
223 & 'WANG= ',f10.6,' (bending)'/
224 & 'WSCLOC= ',f10.6,' (SC local)'/
225 & 'WTOR= ',f10.6,' (torsional)'/
226 & 'WTORD= ',f10.6,' (double torsional)'/
227 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
228 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
229 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
230 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
231 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
232 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
233 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
234 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
235 & 'WSCCOR= ',f10.6,' (SC-backbone torsional correalations)'/
236 & 'WDFAD= ',f10.6,' (DFA distance)'/
237 & 'WDFAT= ',f10.6,' (DFA torsional)'/
238 & 'WDFAN= ',f10.6,' (DFA neighbors)'/
239 & 'WDFAB= ',f10.6,' (DFA beta)'/)
240 if (wcorr4.gt.0.0d0) then
241 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
242 & 'between contact pairs of peptide groups'
243 write (iout,'(2(a,f5.3/))')
244 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
245 & 'Range of quenching the correlation terms:',2*delt_corr
246 else if (wcorr.gt.0.0d0) then
247 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
248 & 'between contact pairs of peptide groups'
250 write (iout,'(a,f8.3)')
251 & 'Scaling factor of 1,4 SC-p interactions:',scal14
252 write (iout,'(a,f8.3)')
253 & 'General scaling factor of SC-p interactions:',scalscp
254 r0_corr=cutoff_corr-delt_corr
256 aad(i,1)=scalscp*aad(i,1)
257 aad(i,2)=scalscp*aad(i,2)
258 bad(i,1)=scalscp*bad(i,1)
259 bad(i,2)=scalscp*bad(i,2)
267 print *,'indpdb=',indpdb,' pdbref=',pdbref
269 C Read sequence if not taken from the pdb file.
270 if (iscode.gt.0) then
271 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
273 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
275 C Convert sequence to numeric code
277 itype(i)=rescode(i,sequence(i),iscode)
279 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
281 & "Glycine is the first full residue, initial dummy deleted"
287 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
289 & "Glycine is the last full residue, terminal dummy deleted"
293 print '(20i4)',(itype(i),i=1,nres)
297 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
299 if (itype(i).eq.21) then
303 else if (itype(i+1).ne.20) then
305 else if (itype(i).ne.20) then
312 write (iout,*) "ITEL"
314 write (iout,*) i,itype(i),itel(i)
317 print *,'Call Read_Bridge.'
319 if (with_dihed_constr) then
321 read (inp,*) ndih_constr
322 if (ndih_constr.gt.0) then
324 write (iout,*) 'FTORS',ftors
325 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
327 & 'There are',ndih_constr,' constraints on phi angles.'
329 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
332 phi0(i)=deg2rad*phi0(i)
333 drange(i)=deg2rad*drange(i)
341 print *,'NNT=',NNT,' NCT=',NCT
342 if (itype(1).eq.21) nnt=2
343 if (itype(nres).eq.21) nct=nct-1
344 if (nstart.lt.nnt) nstart=nnt
345 if (nend.gt.nct .or. nend.eq.0) nend=nct
346 write (iout,*) "nstart",nstart," nend",nend
349 C Juyong:READ init_vars
350 C Initialize variables!
351 C Juyong:READ read_info
352 C READ fragment information!!
353 C both routines should be in dfa.F file!!
355 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
356 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
358 write (iout,*) "Calling init_dfa_vars"
367 write (iout,*) 'init_dfa_vars finished!'
376 write (iout,*) 'read_dfa_info finished!'
385 if (constr_homology.gt.0) then
386 call read_constr_homology(print_homology_restraints)
390 c read(inp,'(a)') pdbfile
391 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
392 c open(ipdbin,file=pdbfile,status='old',err=33)
394 c 33 write (iout,'(a)') 'Error opening PDB file.'
397 c print *,'Begin reading pdb data'
399 c print *,'Finished reading pdb data'
400 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
402 c itype_pdb(i)=itype(i)
405 c write (iout,'(a,i3)') 'nsup=',nsup
407 c if (nsup.le.(nct-nnt+1)) then
408 c do i=0,nct-nnt+1-nsup
409 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
415 c & 'Error - sequences to be superposed do not match.'
418 c do i=0,nsup-(nct-nnt+1)
419 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
421 c nstart_sup=nstart_sup+i
427 c & 'Error - sequences to be superposed do not match.'
430 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
431 c & ' nstart_seq=',nstart_seq
435 write (iout,*) "molread: REFSTR",refstr
437 if (.not.pdbref) then
438 call read_angles(inp,*38)
440 38 write (iout,'(a)') 'Error reading reference structure.'
442 call mp_stopall(Error_Msg)
444 stop 'Error reading reference structure'
456 call contact(print_contact_map,ncont_ref,icont_ref)
458 c Read distance restraints
459 if (constr_dist.gt.0) then
460 call read_dist_constr
465 c-----------------------------------------------------------------------------
466 logical function seq_comp(itypea,itypeb,length)
468 integer length,itypea(length),itypeb(length)
471 if (itypea(i).ne.itypeb(i)) then
479 c-----------------------------------------------------------------------------
480 subroutine read_bridge
481 C Read information about disulfide bridges.
484 include 'COMMON.IOUNITS'
487 include 'COMMON.INTERACT'
488 include 'COMMON.LOCAL'
489 include 'COMMON.NAMES'
490 include 'COMMON.CHAIN'
491 include 'COMMON.FFIELD'
492 include 'COMMON.SBRIDGE'
493 include 'COMMON.HEADER'
494 include 'COMMON.CONTROL'
495 include 'COMMON.TIME1'
497 include 'COMMON.INFO'
500 C Read bridging residues.
501 read (inp,*) ns,(iss(i),i=1,ns)
502 write(iout,*)'ns=',ns
503 C Check whether the specified bridging residues are cystines.
505 if (itype(iss(i)).ne.1) then
506 write (iout,'(2a,i3,a)')
507 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
508 & ' can form a disulfide bridge?!!!'
509 write (*,'(2a,i3,a)')
510 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
511 & ' can form a disulfide bridge?!!!'
513 call mp_stopall(error_msg)
519 C Read preformed bridges.
521 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
524 C Check if the residues involved in bridges are in the specified list of
528 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
529 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
530 write (iout,'(a,i3,a)') 'Disulfide pair',i,
531 & ' contains residues present in other pairs.'
532 write (*,'(a,i3,a)') 'Disulfide pair',i,
533 & ' contains residues present in other pairs.'
535 call mp_stopall(error_msg)
542 if (ihpb(i).eq.iss(j)) goto 10
544 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
547 if (jhpb(i).eq.iss(j)) goto 20
549 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
560 if (ns.gt.0.and.dyn_ss) then
564 forcon(i-nss)=forcon(i)
571 dyn_ss_mask(iss(i))=.true.
572 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
575 print *, "Leaving brigde read"
578 c----------------------------------------------------------------------------
579 subroutine read_angles(kanal,*)
584 include 'COMMON.CHAIN'
585 include 'COMMON.IOUNITS'
587 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
588 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
589 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
590 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
592 theta(i)=deg2rad*theta(i)
593 phi(i)=deg2rad*phi(i)
594 alph(i)=deg2rad*alph(i)
595 omeg(i)=deg2rad*omeg(i)
600 c----------------------------------------------------------------------------
601 subroutine reada(rekord,lancuch,wartosc,default)
603 character*(*) rekord,lancuch
604 double precision wartosc,default
607 iread=index(rekord,lancuch)
612 iread=iread+ilen(lancuch)+1
613 read (rekord(iread:),*) wartosc
616 c----------------------------------------------------------------------------
617 subroutine multreada(rekord,lancuch,tablica,dim,default)
620 double precision tablica(dim),default
621 character*(*) rekord,lancuch
627 iread=index(rekord,lancuch)
628 if (iread.eq.0) return
629 iread=iread+ilen(lancuch)+1
630 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
633 c----------------------------------------------------------------------------
634 subroutine readi(rekord,lancuch,wartosc,default)
636 character*(*) rekord,lancuch
637 integer wartosc,default
640 iread=index(rekord,lancuch)
645 iread=iread+ilen(lancuch)+1
646 read (rekord(iread:),*) wartosc
649 c----------------------------------------------------------------------------
650 subroutine multreadi(rekord,lancuch,tablica,dim,default)
653 integer tablica(dim),default
654 character*(*) rekord,lancuch
661 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
662 if (iread.eq.0) return
663 iread=iread+ilen(lancuch)+1
664 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
667 c----------------------------------------------------------------------------
668 subroutine card_concat(card)
670 include 'COMMON.IOUNITS'
672 character*80 karta,ucase
674 read (inp,'(a)') karta
677 do while (karta(80:80).eq.'&')
678 card=card(:ilen(card)+1)//karta(:79)
679 read (inp,'(a)') karta
682 card=card(:ilen(card)+1)//karta
685 c----------------------------------------------------------------------------
694 include 'COMMON.IOUNITS'
695 include 'COMMON.CONTROL'
696 integer lenpre,lenpot,ilen
698 character*16 cformat,cprint
700 integer lenint,lenout
701 call getenv('INPUT',prefix)
702 call getenv('OUTPUT',prefout)
703 call getenv('INTIN',prefintin)
704 call getenv('COORD',cformat)
705 call getenv('PRINTCOOR',cprint)
706 call getenv('SCRATCHDIR',scratchdir)
709 if (index(ucase(cformat),'CX').gt.0) then
716 lenint=ilen(prefintin)
717 C Get the names and open the input files
718 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
720 write (liczba,'(bz,i3.3)') me
721 outname=prefout(:lenout)//'_clust.out_'//liczba
723 outname=prefout(:lenout)//'_clust.out'
726 intinname=prefintin(:lenint)//'.bx'
727 else if (from_cx) then
728 intinname=prefintin(:lenint)//'.cx'
730 intinname=prefintin(:lenint)//'.int'
732 rmsname=prefintin(:lenint)//'.rms'
733 open (jplot,file=prefout(:ilen(prefout))//'.tex',
735 open (jrms,file=rmsname,status='unknown')
736 open(iout,file=outname,status='unknown')
737 C Get parameter filenames and open the parameter files.
738 call getenv('BONDPAR',bondname)
739 open (ibond,file=bondname,status='old')
740 call getenv('THETPAR',thetname)
741 open (ithep,file=thetname,status='old')
742 call getenv('ROTPAR',rotname)
743 open (irotam,file=rotname,status='old')
744 call getenv('TORPAR',torname)
745 open (itorp,file=torname,status='old')
746 call getenv('TORDPAR',tordname)
747 open (itordp,file=tordname,status='old')
748 call getenv('FOURIER',fouriername)
749 open (ifourier,file=fouriername,status='old')
750 call getenv('ELEPAR',elename)
751 open (ielep,file=elename,status='old')
752 call getenv('SIDEPAR',sidename)
753 open (isidep,file=sidename,status='old')
754 call getenv('SIDEP',sidepname)
755 open (isidep1,file=sidepname,status="old")
756 call getenv('SCCORPAR',sccorname)
757 open (isccor,file=sccorname,status="old")
760 C 8/9/01 In the newest version SCp interaction constants are read from a file
761 C Use -DOLDSCP to use hard-coded constants instead.
763 call getenv('SCPPAR',scpname)
764 open (iscpp,file=scpname,status='old')
768 c-------------------------------------------------------------------------------
769 subroutine read_dist_constr
770 implicit real*8 (a-h,o-z)
775 include 'COMMON.CONTROL'
776 include 'COMMON.CHAIN'
777 include 'COMMON.IOUNITS'
778 include 'COMMON.SBRIDGE'
779 integer ifrag_(2,100),ipair_(2,100)
780 double precision wfrag_(100),wpair_(100)
781 character*500 controlcard
782 logical lprn /.true./
783 logical normalize,next
785 double precision xlink(4,0:4) /
787 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
788 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
789 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
790 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
791 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
792 write (iout,*) "Calling read_dist_constr"
793 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
799 call card_concat(controlcard)
800 next = index(controlcard,"NEXT").gt.0
801 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
802 write (iout,*) "restr_type",restr_type
803 call readi(controlcard,"NFRAG",nfrag_,0)
804 call readi(controlcard,"NFRAG",nfrag_,0)
805 call readi(controlcard,"NPAIR",npair_,0)
806 call readi(controlcard,"NDIST",ndist_,0)
807 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
808 if (restr_type.eq.10)
809 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
810 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
811 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
812 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
813 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
814 normalize = index(controlcard,"NORMALIZE").gt.0
815 write (iout,*) "WBOLTZD",wboltzd
816 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
817 c write (iout,*) "IFRAG"
819 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
821 c write (iout,*) "IPAIR"
823 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
825 if (nfrag_.gt.0) then
827 read(inp,'(a)') pdbfile
829 & "Distance restraints will be constructed from structure ",pdbfile
830 open(ipdbin,file=pdbfile,status='old',err=11)
836 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
837 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
838 & ifrag_(2,i)=nstart_sup+nsup-1
839 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
841 if (wfrag_(i).eq.0.0d0) cycle
842 do j=ifrag_(1,i),ifrag_(2,i)-1
844 c write (iout,*) "j",j," k",k
846 if (restr_type.eq.1) then
852 forcon(nhpb)=wfrag_(i)
853 else if (constr_dist.eq.2) then
854 if (ddjk.le.dist_cut) then
860 forcon(nhpb)=wfrag_(i)
862 else if (restr_type.eq.3) then
868 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
870 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
871 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
876 if (wpair_(i).eq.0.0d0) cycle
884 do j=ifrag_(1,ii),ifrag_(2,ii)
885 do k=ifrag_(1,jj),ifrag_(2,jj)
886 if (restr_type.eq.1) then
892 forcon(nhpb)=wfrag_(i)
893 else if (constr_dist.eq.2) then
894 if (ddjk.le.dist_cut) then
900 forcon(nhpb)=wfrag_(i)
902 else if (restr_type.eq.3) then
908 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
910 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
911 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
917 write (iout,*) "Distance restraints as read from input"
919 if (restr_type.eq.11) then
920 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
921 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
922 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
923 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
926 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
927 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
928 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
929 if (ibecarb(nhpb).gt.0) then
930 ihpb(nhpb)=ihpb(nhpb)+nres
931 jhpb(nhpb)=jhpb(nhpb)+nres
933 else if (constr_dist.eq.10) then
934 c Cross-lonk Markov-like potential
935 call card_concat(controlcard)
936 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
937 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
939 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
940 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
941 if (index(controlcard,"ZL").gt.0) then
943 else if (index(controlcard,"ADH").gt.0) then
945 else if (index(controlcard,"PDH").gt.0) then
947 else if (index(controlcard,"DSS").gt.0) then
952 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
953 & xlink(1,link_type))
954 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
955 & xlink(2,link_type))
956 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
957 & xlink(3,link_type))
958 call reada(controlcard,"SIGMA",forcon(nhpb+1),
959 & xlink(4,link_type))
960 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
961 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
962 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
963 if (forcon(nhpb+1).le.0.0d0 .or.
964 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
967 if (ibecarb(nhpb).gt.0) then
968 ihpb(nhpb)=ihpb(nhpb)+nres
969 jhpb(nhpb)=jhpb(nhpb)+nres
971 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
972 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
973 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
977 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
978 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
979 if (forcon(nhpb+1).gt.0.0d0) then
981 if (dhpb1(nhpb).eq.0.0d0) then
986 if (ibecarb(nhpb).gt.0) then
987 ihpb(nhpb)=ihpb(nhpb)+nres
988 jhpb(nhpb)=jhpb(nhpb)+nres
990 if (dhpb(nhpb).eq.0.0d0)
991 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
993 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
994 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
996 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
997 C if (forcon(nhpb+1).gt.0.0d0) then
999 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1007 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1008 & fordepthmax=fordepth(i)
1011 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1014 if (nhpb.gt.nss) then
1015 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1016 & "The following",nhpb-nss,
1017 & " distance restraints have been imposed:",
1018 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1021 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1022 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1029 11 write (iout,*)"read_dist_restr: error reading reference structure"
1032 c====-------------------------------------------------------------------
1033 subroutine read_constr_homology(lprn)
1035 include 'DIMENSIONS'
1039 include 'COMMON.SETUP'
1040 include 'COMMON.CONTROL'
1041 include 'COMMON.CHAIN'
1042 include 'COMMON.IOUNITS'
1043 include 'COMMON.GEO'
1044 include 'COMMON.INTERACT'
1045 include 'COMMON.NAMES'
1046 include 'COMMON.HOMRESTR'
1048 c For new homol impl
1050 include 'COMMON.VAR'
1051 c include 'include_unres/COMMON.VAR'
1054 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1056 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1057 c & sigma_odl_temp(maxres,maxres,max_template)
1059 character*24 model_ki_dist, model_ki_angle
1060 character*500 controlcard
1061 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1062 integer idomain(max_template,maxres)
1066 logical unres_pdb,liiflag
1068 c FP - Nov. 2014 Temporary specifications for new vars
1070 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1071 double precision, dimension (max_template,maxres) :: rescore
1072 double precision, dimension (max_template,maxres) :: rescore2
1073 character*24 tpl_k_rescore
1074 c -----------------------------------------------------------------
1075 c Reading multiple PDB ref structures and calculation of retraints
1076 c not using pre-computed ones stored in files model_ki_{dist,angle}
1078 c -----------------------------------------------------------------
1081 c Alternative: reading from input
1083 write (iout,*) "BEGIN READ HOMOLOGY INFO"
1090 call card_concat(controlcard)
1091 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1092 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1093 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1094 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1095 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1096 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1097 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1098 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1099 if (homol_nset.gt.1)then
1100 call readi(controlcard,"ISET",iset,1)
1101 call card_concat(controlcard)
1102 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1105 waga_homology(1)=1.0
1109 write(iout,*) "read_constr_homology iset",iset
1110 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1117 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 if (read_homol_frag) then
1149 call read_klapaucjusz
1151 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1152 do k=1,constr_homology
1154 read(inp,'(a)') pdbfile
1155 c write (iout,*) "k ",k," pdbfile ",pdbfile
1156 c Next stament causes error upon compilation (?)
1157 c if(me.eq.king.or. .not. out1file)
1158 c write (iout,'(2a)') 'PDB data will be read from file ',
1159 c & pdbfile(:ilen(pdbfile))
1160 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1161 & pdbfile(:ilen(pdbfile))
1162 open(ipdbin,file=pdbfile,status='old',err=33)
1164 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1165 & pdbfile(:ilen(pdbfile))
1168 c print *,'Begin reading pdb data'
1170 c Files containing res sim or local scores (former containing sigmas)
1173 write(kic2,'(bz,i2.2)') k
1175 tpl_k_rescore="template"//kic2//".sco"
1181 crefjlee(j,i)=c(j,i)
1186 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1187 & (crefjlee(j,i+nres),j=1,3)
1189 write (iout,*) "READ HOMOLOGY INFO"
1190 write (iout,*) "read_constr_homology x: after reading pdb file"
1191 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1192 write (iout,*) "waga_dist",waga_dist
1193 write (iout,*) "waga_angle",waga_angle
1194 write (iout,*) "waga_theta",waga_theta
1195 write (iout,*) "waga_d",waga_d
1196 write (iout,*) "dist_cut",dist_cut
1205 c Distance restraints
1208 C Copy the coordinates from reference coordinates (?)
1212 c write (iout,*) "c(",j,i,") =",c(j,i)
1216 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1218 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1219 open (ientin,file=tpl_k_rescore,status='old')
1220 if (nnt.gt.1) rescore(k,1)=0.0d0
1221 do irec=nnt,nct ! loop for reading res sim
1222 if (read2sigma) then
1223 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1226 idomain(k,i_tmp)=idomain_tmp
1227 rescore(k,i_tmp)=rescore_tmp
1228 rescore2(k,i_tmp)=rescore2_tmp
1229 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
1230 & i_tmp,rescore2_tmp,rescore_tmp,
1234 read (ientin,*,end=1401) rescore_tmp
1236 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1237 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1238 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1243 if (waga_dist.ne.0.0d0) then
1251 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1252 c write (iout,*) k,i,j,distal,dist2_cut
1254 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1255 & .and. distal.le.dist2_cut ) then
1261 c write (iout,*) "k",k
1262 c write (iout,*) "i",i," j",j," constr_homology",
1267 if (read2sigma) then
1270 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1272 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1273 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1274 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1276 if (odl(k,ii).le.dist_cut) then
1277 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1280 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1281 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1283 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1284 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1288 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1291 l_homo(k,ii)=.false.
1298 c Theta, dihedral and SC retraints
1300 if (waga_angle.gt.0.0d0) then
1301 c open (ientin,file=tpl_k_sigma_dih,status='old')
1302 c do irec=1,maxres-3 ! loop for reading sigma_dih
1303 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1304 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1305 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1306 c & sigma_dih(k,i+nnt-1)
1311 if (idomain(k,i).eq.0) then
1315 dih(k,i)=phiref(i) ! right?
1316 c read (ientin,*) sigma_dih(k,i) ! original variant
1317 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1318 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1319 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1320 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1321 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1323 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1324 & rescore(k,i-2)+rescore(k,i-3))/4.0
1325 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1326 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1327 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1328 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1329 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1330 c Instead of res sim other local measure of b/b str reliability possible
1331 if (sigma_dih(k,i).ne.0)
1332 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1333 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1338 if (waga_theta.gt.0.0d0) then
1339 c open (ientin,file=tpl_k_sigma_theta,status='old')
1340 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1341 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1342 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1343 c & sigma_theta(k,i+nnt-1)
1348 do i = nnt+2,nct ! right? without parallel.
1349 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1350 c do i=ithet_start,ithet_end ! with FG parallel.
1351 if (idomain(k,i).eq.0) then
1352 sigma_theta(k,i)=0.0
1355 thetatpl(k,i)=thetaref(i)
1356 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1357 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1358 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1359 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1360 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1361 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1362 & rescore(k,i-2))/3.0
1363 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1364 if (sigma_theta(k,i).ne.0)
1365 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1367 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1368 c rescore(k,i-2) ! right expression ?
1369 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1373 if (waga_d.gt.0.0d0) then
1374 c open (ientin,file=tpl_k_sigma_d,status='old')
1375 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1376 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1377 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1378 c & sigma_d(k,i+nnt-1)
1382 do i = nnt,nct ! right? without parallel.
1383 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1384 c do i=loc_start,loc_end ! with FG parallel.
1385 if (itype(i).eq.10) cycle
1386 if (idomain(k,i).eq.0 ) then
1393 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1394 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1395 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1396 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1397 sigma_d(k,i)=rescore(k,i) ! right expression ?
1398 if (sigma_d(k,i).ne.0)
1399 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1401 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1402 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1403 c read (ientin,*) sigma_d(k,i) ! 1st variant
1404 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1409 c remove distance restraints not used in any model from the list
1410 c shift data in all arrays
1412 if (waga_dist.ne.0.0d0) then
1418 if (ii_in_use(ii).eq.0.and.liiflag) then
1422 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1423 & .not.liiflag.and.ii.eq.lim_odl) then
1424 if (ii.eq.lim_odl) then
1425 iishift=ii-iistart+1
1430 do ki=iistart,lim_odl-iishift
1431 ires_homo(ki)=ires_homo(ki+iishift)
1432 jres_homo(ki)=jres_homo(ki+iishift)
1433 ii_in_use(ki)=ii_in_use(ki+iishift)
1434 do k=1,constr_homology
1435 odl(k,ki)=odl(k,ki+iishift)
1436 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1437 l_homo(k,ki)=l_homo(k,ki+iishift)
1441 lim_odl=lim_odl-iishift
1447 endif ! .not. klapaucjusz
1449 if (constr_homology.gt.0) call homology_partition
1450 if (constr_homology.gt.0) call init_int_table
1451 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1452 cd & "lim_xx=",lim_xx
1453 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1454 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1458 if (.not.lprn) return
1459 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1460 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1461 write (iout,*) "Distance restraints from templates"
1463 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1464 & ii,ires_homo(ii),jres_homo(ii),
1465 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1466 & ki=1,constr_homology)
1468 write (iout,*) "Dihedral angle restraints from templates"
1470 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1471 & (rad2deg*dih(ki,i),
1472 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1474 write (iout,*) "Virtual-bond angle restraints from templates"
1476 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1477 & (rad2deg*thetatpl(ki,i),
1478 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1480 write (iout,*) "SC restraints from templates"
1482 write(iout,'(i5,100(4f8.2,4x))') i,
1483 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1484 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1487 c -----------------------------------------------------------------
1490 c----------------------------------------------------------------------
1491 subroutine read_klapaucjusz
1493 include 'DIMENSIONS'
1497 include 'COMMON.SETUP'
1498 include 'COMMON.CONTROL'
1499 include 'COMMON.CHAIN'
1500 include 'COMMON.IOUNITS'
1501 include 'COMMON.GEO'
1502 include 'COMMON.INTERACT'
1503 include 'COMMON.NAMES'
1504 include 'COMMON.HOMRESTR'
1505 character*256 fragfile
1506 integer ninclust(maxclust),inclust(max_template,maxclust),
1507 & nresclust(maxclust),iresclust(maxres,maxclust)
1510 character*24 model_ki_dist, model_ki_angle
1511 character*500 controlcard
1512 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1513 integer idomain(max_template,maxres)
1514 logical lprn /.true./
1517 logical unres_pdb,liiflag
1520 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1521 double precision, dimension (max_template,maxres) :: rescore
1522 double precision, dimension (max_template,maxres) :: rescore2
1523 character*24 tpl_k_rescore
1526 c For new homol impl
1528 include 'COMMON.VAR'
1530 double precision chomo(3,maxres2+2,max_template)
1531 call getenv("FRAGFILE",fragfile)
1532 open(ientin,file=fragfile,status="old",err=10)
1533 read(ientin,*) constr_homology,nclust
1539 do k=1,constr_homology
1540 read(ientin,'(a)') pdbfile
1541 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1542 & pdbfile(:ilen(pdbfile))
1543 open(ipdbin,file=pdbfile,status='old',err=33)
1545 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1546 & pdbfile(:ilen(pdbfile))
1563 read(ientin,*) ninclust(i),nresclust(i)
1564 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1565 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1568 c Loop over clusters
1571 do ll = 1,ninclust(l)
1579 idomain(k,iresclust(i,l)+1) = 1
1581 idomain(k,iresclust(i,l)) = 1
1585 c Distance restraints
1588 C Copy the coordinates from reference coordinates (?)
1592 c write (iout,*) "c(",j,i,") =",c(j,i)
1595 call int_from_cart(.true.,.false.)
1596 call sc_loc_geom(.false.)
1598 thetaref(i)=theta(i)
1601 if (waga_dist.ne.0.0d0) then
1609 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1610 c write (iout,*) k,i,j,distal,dist2_cut
1612 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1613 & .and. distal.le.dist2_cut ) then
1619 c write (iout,*) "k",k
1620 c write (iout,*) "i",i," j",j," constr_homology",
1625 if (read2sigma) then
1628 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1630 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1631 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1632 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1634 if (odl(k,ii).le.dist_cut) then
1635 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1638 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1639 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1641 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1642 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1646 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1649 c l_homo(k,ii)=.false.
1656 c Theta, dihedral and SC retraints
1658 if (waga_angle.gt.0.0d0) then
1660 if (idomain(k,i).eq.0) then
1661 c sigma_dih(k,i)=0.0
1665 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1666 & rescore(k,i-2)+rescore(k,i-3))/4.0
1667 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1668 c & " sigma_dihed",sigma_dih(k,i)
1669 if (sigma_dih(k,i).ne.0)
1670 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1675 if (waga_theta.gt.0.0d0) then
1677 if (idomain(k,i).eq.0) then
1678 c sigma_theta(k,i)=0.0
1681 thetatpl(k,i)=thetaref(i)
1682 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1683 & rescore(k,i-2))/3.0
1684 if (sigma_theta(k,i).ne.0)
1685 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1689 if (waga_d.gt.0.0d0) then
1691 if (itype(i).eq.10) cycle
1692 if (idomain(k,i).eq.0 ) then
1699 sigma_d(k,i)=rescore(k,i)
1700 if (sigma_d(k,i).ne.0)
1701 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1702 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1708 c remove distance restraints not used in any model from the list
1709 c shift data in all arrays
1711 if (waga_dist.ne.0.0d0) then
1717 if (ii_in_use(ii).eq.0.and.liiflag) then
1721 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1722 & .not.liiflag.and.ii.eq.lim_odl) then
1723 if (ii.eq.lim_odl) then
1724 iishift=ii-iistart+1
1729 do ki=iistart,lim_odl-iishift
1730 ires_homo(ki)=ires_homo(ki+iishift)
1731 jres_homo(ki)=jres_homo(ki+iishift)
1732 ii_in_use(ki)=ii_in_use(ki+iishift)
1733 do k=1,constr_homology
1734 odl(k,ki)=odl(k,ki+iishift)
1735 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1736 l_homo(k,ki)=l_homo(k,ki+iishift)
1740 lim_odl=lim_odl-iishift
1747 10 stop "Error infragment file"
1749 c----------------------------------------------------------------------