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 call readi(controlcard,'RESCALE',rescale_mode,2)
28 call readi(controlcard,'PDBOUT',outpdb,0)
29 call readi(controlcard,'MOL2OUT',outmol2,0)
30 refstr=(index(controlcard,'REFSTR').gt.0)
31 write (iout,*) "REFSTR",refstr
32 pdbref=(index(controlcard,'PDBREF').gt.0)
33 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
34 iscode=index(controlcard,'ONE_LETTER')
35 tree=(index(controlcard,'MAKE_TREE').gt.0)
36 min_var=(index(controlcard,'MINVAR').gt.0)
37 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
38 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
39 call readi(controlcard,'NCUT',ncut,1)
40 call readi(controlcard,'NSTART',nstart,0)
41 call readi(controlcard,'NEND',nend,0)
42 call reada(controlcard,'ECUT',ecut,10.0d0)
43 call reada(controlcard,'PROB',prob_limit,0.99d0)
44 write (iout,*) "Probability limit",prob_limit
45 lgrp=(index(controlcard,'LGRP').gt.0)
46 caonly=(index(controlcard,'CA_ONLY').gt.0)
47 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
48 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
49 call readi(controlcard,'IOPT',iopt,2)
50 lside = index(controlcard,"SIDE").gt.0
51 efree = index(controlcard,"EFREE").gt.0
52 call readi(controlcard,'NTEMP',nT,1)
53 write (iout,*) "nT",nT
54 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
55 write (iout,*) "nT",nT
56 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
58 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
60 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
61 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
62 lprint_int=index(controlcard,"PRINT_INT") .gt.0
63 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
64 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
65 write (iout,*) "with_dihed_constr ",with_dihed_constr,
66 & " CONSTR_DIST",constr_dist
68 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
69 write (iout,*) "with_homology_constr ",with_dihed_constr,
70 & " CONSTR_HOMOLOGY",constr_homology
76 c--------------------------------------------------------------------------
79 C Read molecular data.
83 include 'COMMON.IOUNITS'
86 include 'COMMON.INTERACT'
87 include 'COMMON.LOCAL'
88 include 'COMMON.NAMES'
89 include 'COMMON.CHAIN'
90 include 'COMMON.FFIELD'
91 include 'COMMON.SBRIDGE'
92 include 'COMMON.HEADER'
93 include 'COMMON.CONTROL'
94 include 'COMMON.CONTACTS'
95 include 'COMMON.TIME1'
96 include 'COMMON.TORCNSTR'
100 character*4 sequence(maxres)
101 character*800 weightcard
103 double precision x(maxvar)
104 integer itype_pdb(maxres)
110 C Read weights of the subsequent energy terms.
111 call card_concat(weightcard)
112 call reada(weightcard,'WSC',wsc,1.0d0)
113 call reada(weightcard,'WLONG',wsc,wsc)
114 call reada(weightcard,'WSCP',wscp,1.0d0)
115 call reada(weightcard,'WELEC',welec,1.0D0)
116 call reada(weightcard,'WVDWPP',wvdwpp,welec)
117 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
118 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
119 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
120 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
121 call reada(weightcard,'WTURN3',wturn3,1.0D0)
122 call reada(weightcard,'WTURN4',wturn4,1.0D0)
123 call reada(weightcard,'WTURN6',wturn6,1.0D0)
124 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
125 call reada(weightcard,'WBOND',wbond,1.0D0)
126 call reada(weightcard,'WTOR',wtor,1.0D0)
127 call reada(weightcard,'WTORD',wtor_d,1.0D0)
128 call reada(weightcard,'WANG',wang,1.0D0)
129 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
130 call reada(weightcard,'SCAL14',scal14,0.4D0)
131 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
132 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
133 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
134 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
135 call reada(weightcard,"D0CM",d0cm,3.78d0)
136 call reada(weightcard,"AKCM",akcm,15.1d0)
137 call reada(weightcard,"AKTH",akth,11.0d0)
138 call reada(weightcard,"AKCT",akct,12.0d0)
139 call reada(weightcard,"V1SS",v1ss,-1.08d0)
140 call reada(weightcard,"V2SS",v2ss,7.61d0)
141 call reada(weightcard,"V3SS",v3ss,13.7d0)
142 call reada(weightcard,"EBR",ebr,-5.50D0)
143 if (index(weightcard,'SOFT').gt.0) ipot=6
144 C 12/1/95 Added weight for the multi-body term WCORR
145 call reada(weightcard,'WCORRH',wcorr,1.0D0)
148 dyn_ssbond_ij(i,j)=1.0d300
151 call reada(weightcard,"HT",Ht,0.0D0)
153 ss_depth=ebr/wsc-0.25*eps(1,1)
154 Ht=Ht/wsc-0.25*eps(1,1)
155 akcm=akcm*wstrain/wsc
156 akth=akth*wstrain/wsc
157 akct=akct*wstrain/wsc
158 v1ss=v1ss*wstrain/wsc
159 v2ss=v2ss*wstrain/wsc
160 v3ss=v3ss*wstrain/wsc
162 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
164 write (iout,'(/a)') "Disulfide bridge parameters:"
165 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
166 write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
167 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
168 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
169 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
171 write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
172 if (wcorr4.gt.0.0d0) wcorr=wcorr4
191 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
192 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
193 & wturn4,wturn6,wsccor
194 10 format (/'Energy-term weights (unscaled):'//
195 & 'WSCC= ',f10.6,' (SC-SC)'/
196 & 'WSCP= ',f10.6,' (SC-p)'/
197 & 'WELEC= ',f10.6,' (p-p electr)'/
198 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
199 & 'WBOND= ',f10.6,' (stretching)'/
200 & 'WANG= ',f10.6,' (bending)'/
201 & 'WSCLOC= ',f10.6,' (SC local)'/
202 & 'WTOR= ',f10.6,' (torsional)'/
203 & 'WTORD= ',f10.6,' (double torsional)'/
204 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
205 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
206 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
207 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
208 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
209 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
210 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
211 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
212 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
213 if (wcorr4.gt.0.0d0) then
214 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
215 & 'between contact pairs of peptide groups'
216 write (iout,'(2(a,f5.3/))')
217 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
218 & 'Range of quenching the correlation terms:',2*delt_corr
219 else if (wcorr.gt.0.0d0) then
220 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
221 & 'between contact pairs of peptide groups'
223 write (iout,'(a,f8.3)')
224 & 'Scaling factor of 1,4 SC-p interactions:',scal14
225 write (iout,'(a,f8.3)')
226 & 'General scaling factor of SC-p interactions:',scalscp
227 r0_corr=cutoff_corr-delt_corr
229 aad(i,1)=scalscp*aad(i,1)
230 aad(i,2)=scalscp*aad(i,2)
231 bad(i,1)=scalscp*bad(i,1)
232 bad(i,2)=scalscp*bad(i,2)
236 print *,'indpdb=',indpdb,' pdbref=',pdbref
238 C Read sequence if not taken from the pdb file.
239 if (iscode.gt.0) then
240 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
242 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
244 C Convert sequence to numeric code
246 itype(i)=rescode(i,sequence(i),iscode)
249 print '(20i4)',(itype(i),i=1,nres)
253 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
255 if (itype(i).eq.21) then
259 else if (itype(i+1).ne.20) then
261 else if (itype(i).ne.20) then
268 write (iout,*) "ITEL"
270 write (iout,*) i,itype(i),itel(i)
273 print *,'Call Read_Bridge.'
275 if (with_dihed_constr) then
277 read (inp,*) ndih_constr
278 if (ndih_constr.gt.0) then
280 write (iout,*) 'FTORS',ftors
281 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
283 & 'There are',ndih_constr,' constraints on phi angles.'
285 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
288 phi0(i)=deg2rad*phi0(i)
289 drange(i)=deg2rad*drange(i)
297 print *,'NNT=',NNT,' NCT=',NCT
298 if (itype(1).eq.21) nnt=2
299 if (itype(nres).eq.21) nct=nct-1
300 if (nstart.lt.nnt) nstart=nnt
301 if (nend.gt.nct .or. nend.eq.0) nend=nct
302 write (iout,*) "nstart",nstart," nend",nend
305 if (constr_homology.gt.0) then
306 call read_constr_homology
310 c read(inp,'(a)') pdbfile
311 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
312 c open(ipdbin,file=pdbfile,status='old',err=33)
314 c 33 write (iout,'(a)') 'Error opening PDB file.'
317 c print *,'Begin reading pdb data'
319 c print *,'Finished reading pdb data'
320 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
322 c itype_pdb(i)=itype(i)
325 c write (iout,'(a,i3)') 'nsup=',nsup
327 c if (nsup.le.(nct-nnt+1)) then
328 c do i=0,nct-nnt+1-nsup
329 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
335 c & 'Error - sequences to be superposed do not match.'
338 c do i=0,nsup-(nct-nnt+1)
339 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
341 c nstart_sup=nstart_sup+i
347 c & 'Error - sequences to be superposed do not match.'
350 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
351 c & ' nstart_seq=',nstart_seq
355 write (iout,*) "molread: REFSTR",refstr
357 if (.not.pdbref) then
358 call read_angles(inp,*38)
360 38 write (iout,'(a)') 'Error reading reference structure.'
362 call mp_stopall(Error_Msg)
364 stop 'Error reading reference structure'
376 call contact(.true.,ncont_ref,icont_ref)
378 c Read distance restraints
379 if (constr_dist.gt.0) then
380 call read_dist_constr
385 c-----------------------------------------------------------------------------
386 logical function seq_comp(itypea,itypeb,length)
388 integer length,itypea(length),itypeb(length)
391 if (itypea(i).ne.itypeb(i)) then
399 c-----------------------------------------------------------------------------
400 subroutine read_bridge
401 C Read information about disulfide bridges.
404 include 'COMMON.IOUNITS'
407 include 'COMMON.INTERACT'
408 include 'COMMON.LOCAL'
409 include 'COMMON.NAMES'
410 include 'COMMON.CHAIN'
411 include 'COMMON.FFIELD'
412 include 'COMMON.SBRIDGE'
413 include 'COMMON.HEADER'
414 include 'COMMON.CONTROL'
415 include 'COMMON.TIME1'
417 include 'COMMON.INFO'
420 C Read bridging residues.
421 read (inp,*) ns,(iss(i),i=1,ns)
422 write(iout,*)'ns=',ns
423 C Check whether the specified bridging residues are cystines.
425 if (itype(iss(i)).ne.1) then
426 write (iout,'(2a,i3,a)')
427 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
428 & ' can form a disulfide bridge?!!!'
429 write (*,'(2a,i3,a)')
430 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
431 & ' can form a disulfide bridge?!!!'
433 call mp_stopall(error_msg)
439 C Read preformed bridges.
441 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
444 C Check if the residues involved in bridges are in the specified list of
448 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
449 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
450 write (iout,'(a,i3,a)') 'Disulfide pair',i,
451 & ' contains residues present in other pairs.'
452 write (*,'(a,i3,a)') 'Disulfide pair',i,
453 & ' contains residues present in other pairs.'
455 call mp_stopall(error_msg)
462 if (ihpb(i).eq.iss(j)) goto 10
464 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
467 if (jhpb(i).eq.iss(j)) goto 20
469 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
480 if (ns.gt.0.and.dyn_ss) then
484 forcon(i-nss)=forcon(i)
491 dyn_ss_mask(iss(i))=.true.
492 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
495 print *, "Leaving brigde read"
498 c----------------------------------------------------------------------------
499 subroutine read_angles(kanal,*)
504 include 'COMMON.CHAIN'
505 include 'COMMON.IOUNITS'
507 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
508 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
509 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
510 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
512 theta(i)=deg2rad*theta(i)
513 phi(i)=deg2rad*phi(i)
514 alph(i)=deg2rad*alph(i)
515 omeg(i)=deg2rad*omeg(i)
520 c----------------------------------------------------------------------------
521 subroutine reada(rekord,lancuch,wartosc,default)
523 character*(*) rekord,lancuch
524 double precision wartosc,default
527 iread=index(rekord,lancuch)
532 iread=iread+ilen(lancuch)+1
533 read (rekord(iread:),*) wartosc
536 c----------------------------------------------------------------------------
537 subroutine multreada(rekord,lancuch,tablica,dim,default)
540 double precision tablica(dim),default
541 character*(*) rekord,lancuch
547 iread=index(rekord,lancuch)
548 if (iread.eq.0) return
549 iread=iread+ilen(lancuch)+1
550 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
553 c----------------------------------------------------------------------------
554 subroutine readi(rekord,lancuch,wartosc,default)
556 character*(*) rekord,lancuch
557 integer wartosc,default
560 iread=index(rekord,lancuch)
565 iread=iread+ilen(lancuch)+1
566 read (rekord(iread:),*) wartosc
569 c----------------------------------------------------------------------------
570 subroutine multreadi(rekord,lancuch,tablica,dim,default)
573 integer tablica(dim),default
574 character*(*) rekord,lancuch
581 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
582 if (iread.eq.0) return
583 iread=iread+ilen(lancuch)+1
584 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
587 c----------------------------------------------------------------------------
588 subroutine card_concat(card)
590 include 'COMMON.IOUNITS'
592 character*80 karta,ucase
594 read (inp,'(a)') karta
597 do while (karta(80:80).eq.'&')
598 card=card(:ilen(card)+1)//karta(:79)
599 read (inp,'(a)') karta
602 card=card(:ilen(card)+1)//karta
605 c----------------------------------------------------------------------------
614 include 'COMMON.IOUNITS'
615 include 'COMMON.CONTROL'
616 integer lenpre,lenpot,ilen
618 character*16 cformat,cprint
620 integer lenint,lenout
621 call getenv('INPUT',prefix)
622 call getenv('OUTPUT',prefout)
623 call getenv('INTIN',prefintin)
624 call getenv('COORD',cformat)
625 call getenv('PRINTCOOR',cprint)
626 call getenv('SCRATCHDIR',scratchdir)
629 if (index(ucase(cformat),'CX').gt.0) then
636 lenint=ilen(prefintin)
637 C Get the names and open the input files
638 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
640 write (liczba,'(bz,i3.3)') me
641 outname=prefout(:lenout)//'_clust.out_'//liczba
643 outname=prefout(:lenout)//'_clust.out'
646 intinname=prefintin(:lenint)//'.bx'
647 else if (from_cx) then
648 intinname=prefintin(:lenint)//'.cx'
650 intinname=prefintin(:lenint)//'.int'
652 rmsname=prefintin(:lenint)//'.rms'
653 open (jplot,file=prefout(:ilen(prefout))//'.tex',
655 open (jrms,file=rmsname,status='unknown')
656 open(iout,file=outname,status='unknown')
657 C Get parameter filenames and open the parameter files.
658 call getenv('BONDPAR',bondname)
659 open (ibond,file=bondname,status='old')
660 call getenv('THETPAR',thetname)
661 open (ithep,file=thetname,status='old')
662 call getenv('ROTPAR',rotname)
663 open (irotam,file=rotname,status='old')
664 call getenv('TORPAR',torname)
665 open (itorp,file=torname,status='old')
666 call getenv('TORDPAR',tordname)
667 open (itordp,file=tordname,status='old')
668 call getenv('FOURIER',fouriername)
669 open (ifourier,file=fouriername,status='old')
670 call getenv('ELEPAR',elename)
671 open (ielep,file=elename,status='old')
672 call getenv('SIDEPAR',sidename)
673 open (isidep,file=sidename,status='old')
674 call getenv('SIDEP',sidepname)
675 open (isidep1,file=sidepname,status="old")
676 call getenv('SCCORPAR',sccorname)
677 open (isccor,file=sccorname,status="old")
680 C 8/9/01 In the newest version SCp interaction constants are read from a file
681 C Use -DOLDSCP to use hard-coded constants instead.
683 call getenv('SCPPAR',scpname)
684 open (iscpp,file=scpname,status='old')
688 c-------------------------------------------------------------------------------
689 subroutine read_dist_constr
690 implicit real*8 (a-h,o-z)
692 include 'COMMON.CONTROL'
693 include 'COMMON.CHAIN'
694 include 'COMMON.IOUNITS'
695 include 'COMMON.SBRIDGE'
696 integer ifrag_(2,100),ipair_(2,100)
697 double precision wfrag_(100),wpair_(100)
698 character*500 controlcard
699 c write (iout,*) "Calling read_dist_constr"
700 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
701 call card_concat(controlcard)
703 c write (iout,'(a)') controlcard
704 call readi(controlcard,"NFRAG",nfrag_,0)
705 call readi(controlcard,"NPAIR",npair_,0)
706 call readi(controlcard,"NDIST",ndist_,0)
707 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
708 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
709 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
710 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
711 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
712 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
713 write (iout,*) "IFRAG"
715 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
717 write (iout,*) "IPAIR"
719 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
722 if (.not.refstr .and. nfrag_.gt.0) then
724 & "ERROR: no reference structure to compute distance restraints"
726 & "Restraints must be specified explicitly (NDIST=number)"
729 if (nfrag_.lt.2 .and. npair_.gt.0) then
730 write (iout,*) "ERROR: Less than 2 fragments specified",
731 & " but distance restraints between pairs requested"
736 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
737 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
738 & ifrag_(2,i)=nstart_sup+nsup-1
739 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
741 if (wfrag_(i).gt.0.0d0) then
742 do j=ifrag_(1,i),ifrag_(2,i)-1
744 write (iout,*) "j",j," k",k
746 if (constr_dist.eq.1) then
751 forcon(nhpb)=wfrag_(i)
752 else if (constr_dist.eq.2) then
753 if (ddjk.le.dist_cut) then
758 forcon(nhpb)=wfrag_(i)
765 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
767 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
768 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
774 if (wpair_(i).gt.0.0d0) then
782 do j=ifrag_(1,ii),ifrag_(2,ii)
783 do k=ifrag_(1,jj),ifrag_(2,jj)
787 forcon(nhpb)=wpair_(i)
789 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
790 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
796 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
797 & ibecarb(i),forcon(nhpb+1)
798 if (forcon(nhpb+1).gt.0.0d0) then
800 if (ibecarb(i).gt.0) then
804 if (dhpb(nhpb).eq.0.0d0)
805 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
809 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
810 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
816 c====-------------------------------------------------------------------
817 subroutine read_constr_homology
823 include 'COMMON.CONTROL'
824 include 'COMMON.CHAIN'
825 include 'COMMON.IOUNITS'
826 include 'COMMON.INTERACT'
828 double precision odl_temp,sigma_odl_temp
829 common /przechowalnia/ odl_temp(maxres,maxres,max_template),
830 & sigma_odl_temp(maxres,maxres,max_template)
832 character*24 model_ki_dist, model_ki_angle
833 character*500 controlcard
834 integer ki, i, j, k, l
835 logical lprn /.true./
837 call card_concat(controlcard)
838 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
839 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
840 write (iout,*) "wga_dist",waga_dist," waga_angle",waga_angle
846 do ki=1,constr_homology
847 sigma_odl_temp(i,j,ki)=0.0d0
848 odl_temp(i,j,ki)=0.0d0
853 do ki=1,constr_homology
855 sigma_dih(ki,i)=0.0d0
858 do ki=1,constr_homology
859 write(kic2,'(i2)') ki
860 if (ki.le.9) kic2="0"//kic2(2:2)
862 model_ki_dist="model"//kic2//".dist"
863 model_ki_angle="model"//kic2//".angle"
864 open (ientin,file=model_ki_dist,status='old')
865 do irec=1,maxdim !petla do czytania wiezow na odleglosc
866 read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
867 & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
868 odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
869 sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
870 & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
874 open (ientin,file=model_ki_angle,status='old')
875 do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
876 read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
877 & sigma_dih(ki,i+nnt-1)
878 write (iout,*) i, j, k,l,dih(ki,i+nnt-1),
879 & sigma_dih(ki,i+nnt-1)
880 if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
881 sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
887 write (iout,*) "nnt",nnt," nct",nct
891 c write (iout,*) "i",i," j",j," constr_homology",constr_homology
892 do while (ki.le.constr_homology .and.
893 & sigma_odl_temp(i,j,ki).le.0.0d0)
894 c write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
897 c write (iout,*) "ki",ki
898 if (ki.gt.constr_homology) cycle
902 do ki=1,constr_homology
903 odl(ki,ii)=odl_temp(i,j,ki)
904 sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
909 if (constr_homology.gt.0) call homology_partition
911 if (.not.lprn) return
912 write (iout,*) "Distance restraints from templates"
914 write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
915 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
917 write (iout,*) "Dihedral angle restraints from templates"
919 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
920 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
922 c write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
923 c write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)