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)
144 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
145 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
146 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
147 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
148 if (index(weightcard,'SOFT').gt.0) ipot=6
149 C 12/1/95 Added weight for the multi-body term WCORR
150 call reada(weightcard,'WCORRH',wcorr,1.0D0)
153 dyn_ssbond_ij(i,j)=1.0d300
156 call reada(weightcard,"HT",Ht,0.0D0)
158 ss_depth=ebr/wsc-0.25*eps(1,1)
159 Ht=Ht/wsc-0.25*eps(1,1)
160 akcm=akcm*wstrain/wsc
161 akth=akth*wstrain/wsc
162 akct=akct*wstrain/wsc
163 v1ss=v1ss*wstrain/wsc
164 v2ss=v2ss*wstrain/wsc
165 v3ss=v3ss*wstrain/wsc
167 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
169 write (iout,'(/a)') "Disulfide bridge parameters:"
170 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
171 write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
172 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
173 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
174 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
176 write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
177 if (wcorr4.gt.0.0d0) wcorr=wcorr4
196 weights(22)=wdfa_dist
199 weights(25)=wdfa_beta
200 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
201 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
202 & wturn4,wturn6,wsccor,wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
203 10 format (/'Energy-term weights (unscaled):'//
204 & 'WSCC= ',f10.6,' (SC-SC)'/
205 & 'WSCP= ',f10.6,' (SC-p)'/
206 & 'WELEC= ',f10.6,' (p-p electr)'/
207 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
208 & 'WBOND= ',f10.6,' (stretching)'/
209 & 'WANG= ',f10.6,' (bending)'/
210 & 'WSCLOC= ',f10.6,' (SC local)'/
211 & 'WTOR= ',f10.6,' (torsional)'/
212 & 'WTORD= ',f10.6,' (double torsional)'/
213 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
214 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
215 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
216 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
217 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
218 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
219 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
220 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
221 & 'WSCCOR= ',f10.6,' (SC-backbone torsional correalations)'/
222 & 'WDFAD= ',f10.6,' (DFA distance)'/
223 & 'WDFAT= ',f10.6,' (DFA torsional)'/
224 & 'WDFAN= ',f10.6,' (DFA neighbors)'/
225 & 'WDFAB= ',f10.6,' (DFA beta)'/)
226 if (wcorr4.gt.0.0d0) then
227 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
228 & 'between contact pairs of peptide groups'
229 write (iout,'(2(a,f5.3/))')
230 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
231 & 'Range of quenching the correlation terms:',2*delt_corr
232 else if (wcorr.gt.0.0d0) then
233 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
234 & 'between contact pairs of peptide groups'
236 write (iout,'(a,f8.3)')
237 & 'Scaling factor of 1,4 SC-p interactions:',scal14
238 write (iout,'(a,f8.3)')
239 & 'General scaling factor of SC-p interactions:',scalscp
240 r0_corr=cutoff_corr-delt_corr
242 aad(i,1)=scalscp*aad(i,1)
243 aad(i,2)=scalscp*aad(i,2)
244 bad(i,1)=scalscp*bad(i,1)
245 bad(i,2)=scalscp*bad(i,2)
249 print *,'indpdb=',indpdb,' pdbref=',pdbref
251 C Read sequence if not taken from the pdb file.
252 if (iscode.gt.0) then
253 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
255 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
257 C Convert sequence to numeric code
259 itype(i)=rescode(i,sequence(i),iscode)
262 print '(20i4)',(itype(i),i=1,nres)
266 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
268 if (itype(i).eq.21) then
272 else if (itype(i+1).ne.20) then
274 else if (itype(i).ne.20) then
281 write (iout,*) "ITEL"
283 write (iout,*) i,itype(i),itel(i)
286 print *,'Call Read_Bridge.'
288 if (with_dihed_constr) then
290 read (inp,*) ndih_constr
291 if (ndih_constr.gt.0) then
293 write (iout,*) 'FTORS',ftors
294 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
296 & 'There are',ndih_constr,' constraints on phi angles.'
298 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
301 phi0(i)=deg2rad*phi0(i)
302 drange(i)=deg2rad*drange(i)
310 print *,'NNT=',NNT,' NCT=',NCT
311 if (itype(1).eq.21) nnt=2
312 if (itype(nres).eq.21) nct=nct-1
313 if (nstart.lt.nnt) nstart=nnt
314 if (nend.gt.nct .or. nend.eq.0) nend=nct
315 write (iout,*) "nstart",nstart," nend",nend
318 C Juyong:READ init_vars
319 C Initialize variables!
320 C Juyong:READ read_info
321 C READ fragment information!!
322 C both routines should be in dfa.F file!!
324 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
325 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
326 write (iout,*) "Calling init_dfa_vars"
329 write (iout,*) 'init_dfa_vars finished!'
332 write (iout,*) 'read_dfa_info finished!'
336 if (constr_homology.gt.0) then
337 call read_constr_homology
341 c read(inp,'(a)') pdbfile
342 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
343 c open(ipdbin,file=pdbfile,status='old',err=33)
345 c 33 write (iout,'(a)') 'Error opening PDB file.'
348 c print *,'Begin reading pdb data'
350 c print *,'Finished reading pdb data'
351 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
353 c itype_pdb(i)=itype(i)
356 c write (iout,'(a,i3)') 'nsup=',nsup
358 c if (nsup.le.(nct-nnt+1)) then
359 c do i=0,nct-nnt+1-nsup
360 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
366 c & 'Error - sequences to be superposed do not match.'
369 c do i=0,nsup-(nct-nnt+1)
370 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
372 c nstart_sup=nstart_sup+i
378 c & 'Error - sequences to be superposed do not match.'
381 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
382 c & ' nstart_seq=',nstart_seq
386 write (iout,*) "molread: REFSTR",refstr
388 if (.not.pdbref) then
389 call read_angles(inp,*38)
391 38 write (iout,'(a)') 'Error reading reference structure.'
393 call mp_stopall(Error_Msg)
395 stop 'Error reading reference structure'
407 call contact(.true.,ncont_ref,icont_ref)
409 c Read distance restraints
410 if (constr_dist.gt.0) then
411 call read_dist_constr
416 c-----------------------------------------------------------------------------
417 logical function seq_comp(itypea,itypeb,length)
419 integer length,itypea(length),itypeb(length)
422 if (itypea(i).ne.itypeb(i)) then
430 c-----------------------------------------------------------------------------
431 subroutine read_bridge
432 C Read information about disulfide bridges.
435 include 'COMMON.IOUNITS'
438 include 'COMMON.INTERACT'
439 include 'COMMON.LOCAL'
440 include 'COMMON.NAMES'
441 include 'COMMON.CHAIN'
442 include 'COMMON.FFIELD'
443 include 'COMMON.SBRIDGE'
444 include 'COMMON.HEADER'
445 include 'COMMON.CONTROL'
446 include 'COMMON.TIME1'
448 include 'COMMON.INFO'
451 C Read bridging residues.
452 read (inp,*) ns,(iss(i),i=1,ns)
453 write(iout,*)'ns=',ns
454 C Check whether the specified bridging residues are cystines.
456 if (itype(iss(i)).ne.1) then
457 write (iout,'(2a,i3,a)')
458 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
459 & ' can form a disulfide bridge?!!!'
460 write (*,'(2a,i3,a)')
461 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
462 & ' can form a disulfide bridge?!!!'
464 call mp_stopall(error_msg)
470 C Read preformed bridges.
472 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
475 C Check if the residues involved in bridges are in the specified list of
479 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
480 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
481 write (iout,'(a,i3,a)') 'Disulfide pair',i,
482 & ' contains residues present in other pairs.'
483 write (*,'(a,i3,a)') 'Disulfide pair',i,
484 & ' contains residues present in other pairs.'
486 call mp_stopall(error_msg)
493 if (ihpb(i).eq.iss(j)) goto 10
495 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
498 if (jhpb(i).eq.iss(j)) goto 20
500 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
511 if (ns.gt.0.and.dyn_ss) then
515 forcon(i-nss)=forcon(i)
522 dyn_ss_mask(iss(i))=.true.
523 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
526 print *, "Leaving brigde read"
529 c----------------------------------------------------------------------------
530 subroutine read_angles(kanal,*)
535 include 'COMMON.CHAIN'
536 include 'COMMON.IOUNITS'
538 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
539 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
540 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
541 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
543 theta(i)=deg2rad*theta(i)
544 phi(i)=deg2rad*phi(i)
545 alph(i)=deg2rad*alph(i)
546 omeg(i)=deg2rad*omeg(i)
551 c----------------------------------------------------------------------------
552 subroutine reada(rekord,lancuch,wartosc,default)
554 character*(*) rekord,lancuch
555 double precision wartosc,default
558 iread=index(rekord,lancuch)
563 iread=iread+ilen(lancuch)+1
564 read (rekord(iread:),*) wartosc
567 c----------------------------------------------------------------------------
568 subroutine multreada(rekord,lancuch,tablica,dim,default)
571 double precision tablica(dim),default
572 character*(*) rekord,lancuch
578 iread=index(rekord,lancuch)
579 if (iread.eq.0) return
580 iread=iread+ilen(lancuch)+1
581 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
584 c----------------------------------------------------------------------------
585 subroutine readi(rekord,lancuch,wartosc,default)
587 character*(*) rekord,lancuch
588 integer wartosc,default
591 iread=index(rekord,lancuch)
596 iread=iread+ilen(lancuch)+1
597 read (rekord(iread:),*) wartosc
600 c----------------------------------------------------------------------------
601 subroutine multreadi(rekord,lancuch,tablica,dim,default)
604 integer tablica(dim),default
605 character*(*) rekord,lancuch
612 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
613 if (iread.eq.0) return
614 iread=iread+ilen(lancuch)+1
615 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
618 c----------------------------------------------------------------------------
619 subroutine card_concat(card)
621 include 'COMMON.IOUNITS'
623 character*80 karta,ucase
625 read (inp,'(a)') karta
628 do while (karta(80:80).eq.'&')
629 card=card(:ilen(card)+1)//karta(:79)
630 read (inp,'(a)') karta
633 card=card(:ilen(card)+1)//karta
636 c----------------------------------------------------------------------------
645 include 'COMMON.IOUNITS'
646 include 'COMMON.CONTROL'
647 integer lenpre,lenpot,ilen
649 character*16 cformat,cprint
651 integer lenint,lenout
652 call getenv('INPUT',prefix)
653 call getenv('OUTPUT',prefout)
654 call getenv('INTIN',prefintin)
655 call getenv('COORD',cformat)
656 call getenv('PRINTCOOR',cprint)
657 call getenv('SCRATCHDIR',scratchdir)
660 if (index(ucase(cformat),'CX').gt.0) then
667 lenint=ilen(prefintin)
668 C Get the names and open the input files
669 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
671 write (liczba,'(bz,i3.3)') me
672 outname=prefout(:lenout)//'_clust.out_'//liczba
674 outname=prefout(:lenout)//'_clust.out'
677 intinname=prefintin(:lenint)//'.bx'
678 else if (from_cx) then
679 intinname=prefintin(:lenint)//'.cx'
681 intinname=prefintin(:lenint)//'.int'
683 rmsname=prefintin(:lenint)//'.rms'
684 open (jplot,file=prefout(:ilen(prefout))//'.tex',
686 open (jrms,file=rmsname,status='unknown')
687 open(iout,file=outname,status='unknown')
688 C Get parameter filenames and open the parameter files.
689 call getenv('BONDPAR',bondname)
690 open (ibond,file=bondname,status='old')
691 call getenv('THETPAR',thetname)
692 open (ithep,file=thetname,status='old')
693 call getenv('ROTPAR',rotname)
694 open (irotam,file=rotname,status='old')
695 call getenv('TORPAR',torname)
696 open (itorp,file=torname,status='old')
697 call getenv('TORDPAR',tordname)
698 open (itordp,file=tordname,status='old')
699 call getenv('FOURIER',fouriername)
700 open (ifourier,file=fouriername,status='old')
701 call getenv('ELEPAR',elename)
702 open (ielep,file=elename,status='old')
703 call getenv('SIDEPAR',sidename)
704 open (isidep,file=sidename,status='old')
705 call getenv('SIDEP',sidepname)
706 open (isidep1,file=sidepname,status="old")
707 call getenv('SCCORPAR',sccorname)
708 open (isccor,file=sccorname,status="old")
711 C 8/9/01 In the newest version SCp interaction constants are read from a file
712 C Use -DOLDSCP to use hard-coded constants instead.
714 call getenv('SCPPAR',scpname)
715 open (iscpp,file=scpname,status='old')
719 c-------------------------------------------------------------------------------
720 subroutine read_dist_constr
721 implicit real*8 (a-h,o-z)
723 include 'COMMON.CONTROL'
724 include 'COMMON.CHAIN'
725 include 'COMMON.IOUNITS'
726 include 'COMMON.SBRIDGE'
727 integer ifrag_(2,100),ipair_(2,100)
728 double precision wfrag_(100),wpair_(100)
729 character*500 controlcard
730 c write (iout,*) "Calling read_dist_constr"
731 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
732 call card_concat(controlcard)
734 c write (iout,'(a)') controlcard
735 call readi(controlcard,"NFRAG",nfrag_,0)
736 call readi(controlcard,"NPAIR",npair_,0)
737 call readi(controlcard,"NDIST",ndist_,0)
738 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
739 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
740 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
741 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
742 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
743 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
744 write (iout,*) "IFRAG"
746 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
748 write (iout,*) "IPAIR"
750 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
753 if (.not.refstr .and. nfrag_.gt.0) then
755 & "ERROR: no reference structure to compute distance restraints"
757 & "Restraints must be specified explicitly (NDIST=number)"
760 if (nfrag_.lt.2 .and. npair_.gt.0) then
761 write (iout,*) "ERROR: Less than 2 fragments specified",
762 & " but distance restraints between pairs requested"
767 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
768 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
769 & ifrag_(2,i)=nstart_sup+nsup-1
770 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
772 if (wfrag_(i).gt.0.0d0) then
773 do j=ifrag_(1,i),ifrag_(2,i)-1
775 write (iout,*) "j",j," k",k
777 if (constr_dist.eq.1) then
782 forcon(nhpb)=wfrag_(i)
783 else if (constr_dist.eq.2) then
784 if (ddjk.le.dist_cut) then
789 forcon(nhpb)=wfrag_(i)
796 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
798 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
799 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
805 if (wpair_(i).gt.0.0d0) then
813 do j=ifrag_(1,ii),ifrag_(2,ii)
814 do k=ifrag_(1,jj),ifrag_(2,jj)
818 forcon(nhpb)=wpair_(i)
820 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
821 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
827 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
828 & ibecarb(i),forcon(nhpb+1)
829 if (forcon(nhpb+1).gt.0.0d0) then
831 if (ibecarb(i).gt.0) then
835 if (dhpb(nhpb).eq.0.0d0)
836 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
840 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
841 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
847 c====-------------------------------------------------------------------
848 subroutine read_constr_homology
854 include 'COMMON.CONTROL'
855 include 'COMMON.CHAIN'
856 include 'COMMON.IOUNITS'
857 include 'COMMON.INTERACT'
859 double precision odl_temp,sigma_odl_temp
860 common /przechowalnia/ odl_temp(maxres,maxres,max_template),
861 & sigma_odl_temp(maxres,maxres,max_template)
863 character*24 model_ki_dist, model_ki_angle
864 character*500 controlcard
865 integer ki, i, j, k, l
866 logical lprn /.true./
868 call card_concat(controlcard)
869 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
870 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
871 write (iout,*) "wga_dist",waga_dist," waga_angle",waga_angle
877 do ki=1,constr_homology
878 sigma_odl_temp(i,j,ki)=0.0d0
879 odl_temp(i,j,ki)=0.0d0
884 do ki=1,constr_homology
886 sigma_dih(ki,i)=0.0d0
889 do ki=1,constr_homology
890 write(kic2,'(i2)') ki
891 if (ki.le.9) kic2="0"//kic2(2:2)
893 model_ki_dist="model"//kic2//".dist"
894 model_ki_angle="model"//kic2//".angle"
895 open (imol2,file=model_ki_dist,status='old')
896 do irec=1,maxdim !petla do czytania wiezow na odleglosc
897 read (imol2,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
898 & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
899 odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
900 sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
901 & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
905 open (imol2,file=model_ki_angle,status='old')
906 do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
907 read (imol2,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
908 & sigma_dih(ki,i+nnt-1)
909 write (iout,*) i, j, k,l,dih(ki,i+nnt-1),
910 & sigma_dih(ki,i+nnt-1)
911 if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
912 sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
918 write (iout,*) "nnt",nnt," nct",nct
922 c write (iout,*) "i",i," j",j," constr_homology",constr_homology
923 do while (ki.le.constr_homology .and.
924 & sigma_odl_temp(i,j,ki).le.0.0d0)
925 c write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
928 c write (iout,*) "ki",ki
929 if (ki.gt.constr_homology) cycle
933 do ki=1,constr_homology
934 odl(ki,ii)=odl_temp(i,j,ki)
935 sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
940 if (constr_homology.gt.0) call homology_partition
942 if (.not.lprn) return
943 write (iout,*) "Distance restraints from templates"
945 write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
946 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
948 write (iout,*) "Dihedral angle restraints from templates"
950 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
951 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
953 c write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
954 c write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)