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
88 c--------------------------------------------------------------------------
91 C Read molecular data.
95 include 'COMMON.IOUNITS'
98 include 'COMMON.INTERACT'
99 include 'COMMON.LOCAL'
100 include 'COMMON.NAMES'
101 include 'COMMON.CHAIN'
102 include 'COMMON.FFIELD'
103 include 'COMMON.SBRIDGE'
104 include 'COMMON.HEADER'
105 include 'COMMON.CONTROL'
106 include 'COMMON.CONTACTS'
107 include 'COMMON.TIME1'
108 include 'COMMON.TORCNSTR'
110 include 'COMMON.INFO'
112 character*4 sequence(maxres)
113 character*800 weightcard
115 double precision x(maxvar)
116 integer itype_pdb(maxres)
119 write (iout,*) " MOLREAD: NRES",NRES
123 C Read weights of the subsequent energy terms.
124 call card_concat(weightcard)
125 call reada(weightcard,'WSC',wsc,1.0d0)
126 call reada(weightcard,'WLONG',wsc,wsc)
127 call reada(weightcard,'WSCP',wscp,1.0d0)
128 call reada(weightcard,'WELEC',welec,1.0D0)
129 call reada(weightcard,'WVDWPP',wvdwpp,welec)
130 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
131 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
132 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
133 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
134 call reada(weightcard,'WTURN3',wturn3,1.0D0)
135 call reada(weightcard,'WTURN4',wturn4,1.0D0)
136 call reada(weightcard,'WTURN6',wturn6,1.0D0)
137 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
138 call reada(weightcard,'WBOND',wbond,1.0D0)
139 call reada(weightcard,'WTOR',wtor,1.0D0)
140 call reada(weightcard,'WTORD',wtor_d,1.0D0)
141 call reada(weightcard,'WANG',wang,1.0D0)
142 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
143 call reada(weightcard,'SCAL14',scal14,0.4D0)
144 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
145 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
146 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
147 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
148 call reada(weightcard,"D0CM",d0cm,3.78d0)
149 call reada(weightcard,"AKCM",akcm,15.1d0)
150 call reada(weightcard,"AKTH",akth,11.0d0)
151 call reada(weightcard,"AKCT",akct,12.0d0)
152 call reada(weightcard,"V1SS",v1ss,-1.08d0)
153 call reada(weightcard,"V2SS",v2ss,7.61d0)
154 call reada(weightcard,"V3SS",v3ss,13.7d0)
155 call reada(weightcard,"EBR",ebr,-5.50D0)
157 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
158 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
159 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
160 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
161 if (index(weightcard,'SOFT').gt.0) ipot=6
162 C 12/1/95 Added weight for the multi-body term WCORR
163 call reada(weightcard,'WCORRH',wcorr,1.0D0)
166 dyn_ssbond_ij(i,j)=1.0d300
169 call reada(weightcard,"HT",Ht,0.0D0)
171 ss_depth=ebr/wsc-0.25*eps(1,1)
172 Ht=Ht/wsc-0.25*eps(1,1)
173 akcm=akcm*wstrain/wsc
174 akth=akth*wstrain/wsc
175 akct=akct*wstrain/wsc
176 v1ss=v1ss*wstrain/wsc
177 v2ss=v2ss*wstrain/wsc
178 v3ss=v3ss*wstrain/wsc
180 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
182 write (iout,'(/a)') "Disulfide bridge parameters:"
183 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
184 write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
185 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
186 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
187 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
189 write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
190 if (wcorr4.gt.0.0d0) wcorr=wcorr4
209 weights(22)=wdfa_dist
212 weights(25)=wdfa_beta
213 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
214 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
215 & wturn4,wturn6,wsccor,wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
216 10 format (/'Energy-term weights (unscaled):'//
217 & 'WSCC= ',f10.6,' (SC-SC)'/
218 & 'WSCP= ',f10.6,' (SC-p)'/
219 & 'WELEC= ',f10.6,' (p-p electr)'/
220 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
221 & 'WBOND= ',f10.6,' (stretching)'/
222 & 'WANG= ',f10.6,' (bending)'/
223 & 'WSCLOC= ',f10.6,' (SC local)'/
224 & 'WTOR= ',f10.6,' (torsional)'/
225 & 'WTORD= ',f10.6,' (double torsional)'/
226 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
227 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
228 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
229 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
230 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
231 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
232 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
233 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
234 & 'WSCCOR= ',f10.6,' (SC-backbone torsional correalations)'/
235 & 'WDFAD= ',f10.6,' (DFA distance)'/
236 & 'WDFAT= ',f10.6,' (DFA torsional)'/
237 & 'WDFAN= ',f10.6,' (DFA neighbors)'/
238 & 'WDFAB= ',f10.6,' (DFA beta)'/)
239 if (wcorr4.gt.0.0d0) then
240 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
241 & 'between contact pairs of peptide groups'
242 write (iout,'(2(a,f5.3/))')
243 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
244 & 'Range of quenching the correlation terms:',2*delt_corr
245 else if (wcorr.gt.0.0d0) then
246 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
247 & 'between contact pairs of peptide groups'
249 write (iout,'(a,f8.3)')
250 & 'Scaling factor of 1,4 SC-p interactions:',scal14
251 write (iout,'(a,f8.3)')
252 & 'General scaling factor of SC-p interactions:',scalscp
253 r0_corr=cutoff_corr-delt_corr
255 aad(i,1)=scalscp*aad(i,1)
256 aad(i,2)=scalscp*aad(i,2)
257 bad(i,1)=scalscp*bad(i,1)
258 bad(i,2)=scalscp*bad(i,2)
266 print *,'indpdb=',indpdb,' pdbref=',pdbref
268 C Read sequence if not taken from the pdb file.
269 if (iscode.gt.0) then
270 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
272 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
274 C Convert sequence to numeric code
276 itype(i)=rescode(i,sequence(i),iscode)
278 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
280 & "Glycine is the first full residue, initial dummy deleted"
286 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
288 & "Glycine is the last full residue, terminal dummy deleted"
292 print '(20i4)',(itype(i),i=1,nres)
296 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
298 if (itype(i).eq.21) then
302 else if (itype(i+1).ne.20) then
304 else if (itype(i).ne.20) then
311 write (iout,*) "ITEL"
313 write (iout,*) i,itype(i),itel(i)
316 print *,'Call Read_Bridge.'
318 if (with_dihed_constr) then
320 read (inp,*) ndih_constr
321 if (ndih_constr.gt.0) then
323 write (iout,*) 'FTORS',ftors
324 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
326 & 'There are',ndih_constr,' constraints on phi angles.'
328 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
331 phi0(i)=deg2rad*phi0(i)
332 drange(i)=deg2rad*drange(i)
340 print *,'NNT=',NNT,' NCT=',NCT
341 if (itype(1).eq.21) nnt=2
342 if (itype(nres).eq.21) nct=nct-1
343 if (nstart.lt.nnt) nstart=nnt
344 if (nend.gt.nct .or. nend.eq.0) nend=nct
345 write (iout,*) "nstart",nstart," nend",nend
348 C Juyong:READ init_vars
349 C Initialize variables!
350 C Juyong:READ read_info
351 C READ fragment information!!
352 C both routines should be in dfa.F file!!
354 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
355 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
357 write (iout,*) "Calling init_dfa_vars"
366 write (iout,*) 'init_dfa_vars finished!'
375 write (iout,*) 'read_dfa_info finished!'
384 if (constr_homology.gt.0) then
385 call read_constr_homology(print_homology_restraints)
389 c read(inp,'(a)') pdbfile
390 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
391 c open(ipdbin,file=pdbfile,status='old',err=33)
393 c 33 write (iout,'(a)') 'Error opening PDB file.'
396 c print *,'Begin reading pdb data'
398 c print *,'Finished reading pdb data'
399 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
401 c itype_pdb(i)=itype(i)
404 c write (iout,'(a,i3)') 'nsup=',nsup
406 c if (nsup.le.(nct-nnt+1)) then
407 c do i=0,nct-nnt+1-nsup
408 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
414 c & 'Error - sequences to be superposed do not match.'
417 c do i=0,nsup-(nct-nnt+1)
418 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
420 c nstart_sup=nstart_sup+i
426 c & 'Error - sequences to be superposed do not match.'
429 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
430 c & ' nstart_seq=',nstart_seq
434 write (iout,*) "molread: REFSTR",refstr
436 if (.not.pdbref) then
437 call read_angles(inp,*38)
439 38 write (iout,'(a)') 'Error reading reference structure.'
441 call mp_stopall(Error_Msg)
443 stop 'Error reading reference structure'
455 call contact(print_contact_map,ncont_ref,icont_ref)
457 c Read distance restraints
458 if (constr_dist.gt.0) then
459 call read_dist_constr
464 c-----------------------------------------------------------------------------
465 logical function seq_comp(itypea,itypeb,length)
467 integer length,itypea(length),itypeb(length)
470 if (itypea(i).ne.itypeb(i)) then
478 c-----------------------------------------------------------------------------
479 subroutine read_bridge
480 C Read information about disulfide bridges.
483 include 'COMMON.IOUNITS'
486 include 'COMMON.INTERACT'
487 include 'COMMON.LOCAL'
488 include 'COMMON.NAMES'
489 include 'COMMON.CHAIN'
490 include 'COMMON.FFIELD'
491 include 'COMMON.SBRIDGE'
492 include 'COMMON.HEADER'
493 include 'COMMON.CONTROL'
494 include 'COMMON.TIME1'
496 include 'COMMON.INFO'
499 C Read bridging residues.
500 read (inp,*) ns,(iss(i),i=1,ns)
501 write(iout,*)'ns=',ns
502 C Check whether the specified bridging residues are cystines.
504 if (itype(iss(i)).ne.1) then
505 write (iout,'(2a,i3,a)')
506 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
507 & ' can form a disulfide bridge?!!!'
508 write (*,'(2a,i3,a)')
509 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
510 & ' can form a disulfide bridge?!!!'
512 call mp_stopall(error_msg)
518 C Read preformed bridges.
520 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
523 C Check if the residues involved in bridges are in the specified list of
527 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
528 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
529 write (iout,'(a,i3,a)') 'Disulfide pair',i,
530 & ' contains residues present in other pairs.'
531 write (*,'(a,i3,a)') 'Disulfide pair',i,
532 & ' contains residues present in other pairs.'
534 call mp_stopall(error_msg)
541 if (ihpb(i).eq.iss(j)) goto 10
543 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
546 if (jhpb(i).eq.iss(j)) goto 20
548 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
559 if (ns.gt.0.and.dyn_ss) then
563 forcon(i-nss)=forcon(i)
570 dyn_ss_mask(iss(i))=.true.
571 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
574 print *, "Leaving brigde read"
577 c----------------------------------------------------------------------------
578 subroutine read_angles(kanal,*)
583 include 'COMMON.CHAIN'
584 include 'COMMON.IOUNITS'
586 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
587 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
588 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
589 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
591 theta(i)=deg2rad*theta(i)
592 phi(i)=deg2rad*phi(i)
593 alph(i)=deg2rad*alph(i)
594 omeg(i)=deg2rad*omeg(i)
599 c----------------------------------------------------------------------------
600 subroutine reada(rekord,lancuch,wartosc,default)
602 character*(*) rekord,lancuch
603 double precision wartosc,default
606 iread=index(rekord,lancuch)
611 iread=iread+ilen(lancuch)+1
612 read (rekord(iread:),*) wartosc
615 c----------------------------------------------------------------------------
616 subroutine multreada(rekord,lancuch,tablica,dim,default)
619 double precision tablica(dim),default
620 character*(*) rekord,lancuch
626 iread=index(rekord,lancuch)
627 if (iread.eq.0) return
628 iread=iread+ilen(lancuch)+1
629 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
632 c----------------------------------------------------------------------------
633 subroutine readi(rekord,lancuch,wartosc,default)
635 character*(*) rekord,lancuch
636 integer wartosc,default
639 iread=index(rekord,lancuch)
644 iread=iread+ilen(lancuch)+1
645 read (rekord(iread:),*) wartosc
648 c----------------------------------------------------------------------------
649 subroutine multreadi(rekord,lancuch,tablica,dim,default)
652 integer tablica(dim),default
653 character*(*) rekord,lancuch
660 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
661 if (iread.eq.0) return
662 iread=iread+ilen(lancuch)+1
663 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
666 c----------------------------------------------------------------------------
667 subroutine card_concat(card)
669 include 'COMMON.IOUNITS'
671 character*80 karta,ucase
673 read (inp,'(a)') karta
676 do while (karta(80:80).eq.'&')
677 card=card(:ilen(card)+1)//karta(:79)
678 read (inp,'(a)') karta
681 card=card(:ilen(card)+1)//karta
684 c----------------------------------------------------------------------------
693 include 'COMMON.IOUNITS'
694 include 'COMMON.CONTROL'
695 integer lenpre,lenpot,ilen
697 character*16 cformat,cprint
699 integer lenint,lenout
700 call getenv('INPUT',prefix)
701 call getenv('OUTPUT',prefout)
702 call getenv('INTIN',prefintin)
703 call getenv('COORD',cformat)
704 call getenv('PRINTCOOR',cprint)
705 call getenv('SCRATCHDIR',scratchdir)
708 if (index(ucase(cformat),'CX').gt.0) then
715 lenint=ilen(prefintin)
716 C Get the names and open the input files
717 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
719 write (liczba,'(bz,i3.3)') me
720 outname=prefout(:lenout)//'_clust.out_'//liczba
722 outname=prefout(:lenout)//'_clust.out'
725 intinname=prefintin(:lenint)//'.bx'
726 else if (from_cx) then
727 intinname=prefintin(:lenint)//'.cx'
729 intinname=prefintin(:lenint)//'.int'
731 rmsname=prefintin(:lenint)//'.rms'
732 open (jplot,file=prefout(:ilen(prefout))//'.tex',
734 open (jrms,file=rmsname,status='unknown')
735 open(iout,file=outname,status='unknown')
736 C Get parameter filenames and open the parameter files.
737 call getenv('BONDPAR',bondname)
738 open (ibond,file=bondname,status='old')
739 call getenv('THETPAR',thetname)
740 open (ithep,file=thetname,status='old')
741 call getenv('ROTPAR',rotname)
742 open (irotam,file=rotname,status='old')
743 call getenv('TORPAR',torname)
744 open (itorp,file=torname,status='old')
745 call getenv('TORDPAR',tordname)
746 open (itordp,file=tordname,status='old')
747 call getenv('FOURIER',fouriername)
748 open (ifourier,file=fouriername,status='old')
749 call getenv('ELEPAR',elename)
750 open (ielep,file=elename,status='old')
751 call getenv('SIDEPAR',sidename)
752 open (isidep,file=sidename,status='old')
753 call getenv('SIDEP',sidepname)
754 open (isidep1,file=sidepname,status="old")
755 call getenv('SCCORPAR',sccorname)
756 open (isccor,file=sccorname,status="old")
759 C 8/9/01 In the newest version SCp interaction constants are read from a file
760 C Use -DOLDSCP to use hard-coded constants instead.
762 call getenv('SCPPAR',scpname)
763 open (iscpp,file=scpname,status='old')
767 c-------------------------------------------------------------------------------
768 subroutine read_dist_constr
769 implicit real*8 (a-h,o-z)
771 include 'COMMON.CONTROL'
772 include 'COMMON.CHAIN'
773 include 'COMMON.IOUNITS'
774 include 'COMMON.SBRIDGE'
775 integer ifrag_(2,100),ipair_(2,100)
776 double precision wfrag_(100),wpair_(100)
777 character*500 controlcard
778 c write (iout,*) "Calling read_dist_constr"
779 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
780 call card_concat(controlcard)
782 c write (iout,'(a)') controlcard
783 call readi(controlcard,"NFRAG",nfrag_,0)
784 call readi(controlcard,"NPAIR",npair_,0)
785 call readi(controlcard,"NDIST",ndist_,0)
786 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
787 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
788 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
789 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
790 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
791 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
792 write (iout,*) "IFRAG"
794 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
796 write (iout,*) "IPAIR"
798 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
805 if (.not.refstr .and. nfrag_.gt.0) then
807 & "ERROR: no reference structure to compute distance restraints"
809 & "Restraints must be specified explicitly (NDIST=number)"
812 if (nfrag_.lt.2 .and. npair_.gt.0) then
813 write (iout,*) "ERROR: Less than 2 fragments specified",
814 & " but distance restraints between pairs requested"
823 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
824 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
825 & ifrag_(2,i)=nstart_sup+nsup-1
826 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
828 if (wfrag_(i).gt.0.0d0) then
829 do j=ifrag_(1,i),ifrag_(2,i)-1
831 write (iout,*) "j",j," k",k
833 if (constr_dist.eq.1) then
838 forcon(nhpb)=wfrag_(i)
839 else if (constr_dist.eq.2) then
840 if (ddjk.le.dist_cut) then
845 forcon(nhpb)=wfrag_(i)
852 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
854 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
855 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
861 if (wpair_(i).gt.0.0d0) then
869 do j=ifrag_(1,ii),ifrag_(2,ii)
870 do k=ifrag_(1,jj),ifrag_(2,jj)
874 forcon(nhpb)=wpair_(i)
876 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
877 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
883 if (constr_dist.eq.11) then
884 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
885 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
886 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
888 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
889 & ibecarb(i),forcon(nhpb+1)
891 if (forcon(nhpb+1).gt.0.0d0) then
893 if (ibecarb(i).gt.0) then
897 if (dhpb(nhpb).eq.0.0d0)
898 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
902 if (constr_dist.eq.11) then
903 write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ",
904 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i),
907 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
908 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
919 c====-------------------------------------------------------------------
920 subroutine read_constr_homology(lprn)
926 include 'COMMON.SETUP'
927 include 'COMMON.CONTROL'
928 include 'COMMON.CHAIN'
929 include 'COMMON.IOUNITS'
931 include 'COMMON.INTERACT'
932 include 'COMMON.NAMES'
933 include 'COMMON.HOMRESTR'
938 c include 'include_unres/COMMON.VAR'
941 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
943 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
944 c & sigma_odl_temp(maxres,maxres,max_template)
946 character*24 model_ki_dist, model_ki_angle
947 character*500 controlcard
948 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
949 integer idomain(max_template,maxres)
953 logical unres_pdb,liiflag
955 c FP - Nov. 2014 Temporary specifications for new vars
957 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
958 double precision, dimension (max_template,maxres) :: rescore
959 double precision, dimension (max_template,maxres) :: rescore2
960 character*24 tpl_k_rescore
961 c -----------------------------------------------------------------
962 c Reading multiple PDB ref structures and calculation of retraints
963 c not using pre-computed ones stored in files model_ki_{dist,angle}
965 c -----------------------------------------------------------------
968 c Alternative: reading from input
970 write (iout,*) "BEGIN READ HOMOLOGY INFO"
977 call card_concat(controlcard)
978 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
979 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
980 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
981 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
982 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
983 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
984 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
985 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
986 if (homol_nset.gt.1)then
987 call readi(controlcard,"ISET",iset,1)
988 call card_concat(controlcard)
989 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
996 write(iout,*) "read_constr_homology iset",iset
997 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1004 cd write (iout,*) "nnt",nnt," nct",nct
1011 c Reading HM global scores (prob not required)
1014 do k=1,constr_homology
1018 c open (4,file="HMscore")
1019 c do k=1,constr_homology
1020 c read (4,*,end=521) hmscore_tmp
1021 c hmscore(k)=hmscore_tmp ! Another transformation can be used
1022 c write(*,*) "Model", k, ":", hmscore(k)
1033 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1035 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1036 do k=1,constr_homology
1038 read(inp,'(a)') pdbfile
1039 c write (iout,*) "k ",k," pdbfile ",pdbfile
1040 c Next stament causes error upon compilation (?)
1041 c if(me.eq.king.or. .not. out1file)
1042 c write (iout,'(2a)') 'PDB data will be read from file ',
1043 c & pdbfile(:ilen(pdbfile))
1044 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1045 & pdbfile(:ilen(pdbfile))
1046 open(ipdbin,file=pdbfile,status='old',err=33)
1048 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1049 & pdbfile(:ilen(pdbfile))
1052 c print *,'Begin reading pdb data'
1054 c Files containing res sim or local scores (former containing sigmas)
1057 write(kic2,'(bz,i2.2)') k
1059 tpl_k_rescore="template"//kic2//".sco"
1065 crefjlee(j,i)=c(j,i)
1070 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1071 & (crefjlee(j,i+nres),j=1,3)
1073 write (iout,*) "READ HOMOLOGY INFO"
1074 write (iout,*) "read_constr_homology x: after reading pdb file"
1075 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1076 write (iout,*) "waga_dist",waga_dist
1077 write (iout,*) "waga_angle",waga_angle
1078 write (iout,*) "waga_theta",waga_theta
1079 write (iout,*) "waga_d",waga_d
1080 write (iout,*) "dist_cut",dist_cut
1089 c Distance restraints
1092 C Copy the coordinates from reference coordinates (?)
1096 c write (iout,*) "c(",j,i,") =",c(j,i)
1100 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1102 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1103 open (ientin,file=tpl_k_rescore,status='old')
1104 if (nnt.gt.1) rescore(k,1)=0.0d0
1105 do irec=nnt,nct ! loop for reading res sim
1106 if (read2sigma) then
1107 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1110 idomain(k,i_tmp)=idomain_tmp
1111 rescore(k,i_tmp)=rescore_tmp
1112 rescore2(k,i_tmp)=rescore2_tmp
1113 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
1114 & i_tmp,rescore2_tmp,rescore_tmp,
1118 read (ientin,*,end=1401) rescore_tmp
1120 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1121 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1122 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1127 if (waga_dist.ne.0.0d0) then
1135 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1136 c write (iout,*) k,i,j,distal,dist2_cut
1138 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1139 & .and. distal.le.dist2_cut ) then
1145 c write (iout,*) "k",k
1146 c write (iout,*) "i",i," j",j," constr_homology",
1151 if (read2sigma) then
1154 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1156 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1157 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1158 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1160 if (odl(k,ii).le.dist_cut) then
1161 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1164 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1165 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1167 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1168 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1172 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1175 l_homo(k,ii)=.false.
1182 c Theta, dihedral and SC retraints
1184 if (waga_angle.gt.0.0d0) then
1185 c open (ientin,file=tpl_k_sigma_dih,status='old')
1186 c do irec=1,maxres-3 ! loop for reading sigma_dih
1187 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1188 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1189 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1190 c & sigma_dih(k,i+nnt-1)
1195 if (idomain(k,i).eq.0) then
1199 dih(k,i)=phiref(i) ! right?
1200 c read (ientin,*) sigma_dih(k,i) ! original variant
1201 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1202 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1203 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1204 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1205 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1207 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1208 & rescore(k,i-2)+rescore(k,i-3))/4.0
1209 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1210 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1211 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1212 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1213 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1214 c Instead of res sim other local measure of b/b str reliability possible
1215 if (sigma_dih(k,i).ne.0)
1216 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1217 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1222 if (waga_theta.gt.0.0d0) then
1223 c open (ientin,file=tpl_k_sigma_theta,status='old')
1224 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1225 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1226 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1227 c & sigma_theta(k,i+nnt-1)
1232 do i = nnt+2,nct ! right? without parallel.
1233 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1234 c do i=ithet_start,ithet_end ! with FG parallel.
1235 if (idomain(k,i).eq.0) then
1236 sigma_theta(k,i)=0.0
1239 thetatpl(k,i)=thetaref(i)
1240 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1241 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1242 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1243 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1244 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1245 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1246 & rescore(k,i-2))/3.0
1247 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1248 if (sigma_theta(k,i).ne.0)
1249 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1251 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1252 c rescore(k,i-2) ! right expression ?
1253 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1257 if (waga_d.gt.0.0d0) then
1258 c open (ientin,file=tpl_k_sigma_d,status='old')
1259 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1260 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1261 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1262 c & sigma_d(k,i+nnt-1)
1266 do i = nnt,nct ! right? without parallel.
1267 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1268 c do i=loc_start,loc_end ! with FG parallel.
1269 if (itype(i).eq.10) cycle
1270 if (idomain(k,i).eq.0 ) then
1277 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1278 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1279 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1280 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1281 sigma_d(k,i)=rescore(k,i) ! right expression ?
1282 if (sigma_d(k,i).ne.0)
1283 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1285 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1286 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1287 c read (ientin,*) sigma_d(k,i) ! 1st variant
1288 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1293 c remove distance restraints not used in any model from the list
1294 c shift data in all arrays
1296 if (waga_dist.ne.0.0d0) then
1302 if (ii_in_use(ii).eq.0.and.liiflag) then
1306 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1307 & .not.liiflag.and.ii.eq.lim_odl) then
1308 if (ii.eq.lim_odl) then
1309 iishift=ii-iistart+1
1314 do ki=iistart,lim_odl-iishift
1315 ires_homo(ki)=ires_homo(ki+iishift)
1316 jres_homo(ki)=jres_homo(ki+iishift)
1317 ii_in_use(ki)=ii_in_use(ki+iishift)
1318 do k=1,constr_homology
1319 odl(k,ki)=odl(k,ki+iishift)
1320 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1321 l_homo(k,ki)=l_homo(k,ki+iishift)
1325 lim_odl=lim_odl-iishift
1330 if (constr_homology.gt.0) call homology_partition
1331 if (constr_homology.gt.0) call init_int_table
1332 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1333 cd & "lim_xx=",lim_xx
1334 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1335 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1339 if (.not.lprn) return
1340 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1341 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1342 write (iout,*) "Distance restraints from templates"
1344 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1345 & ii,ires_homo(ii),jres_homo(ii),
1346 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1347 & ki=1,constr_homology)
1349 write (iout,*) "Dihedral angle restraints from templates"
1351 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1352 & (rad2deg*dih(ki,i),
1353 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1355 write (iout,*) "Virtual-bond angle restraints from templates"
1357 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1358 & (rad2deg*thetatpl(ki,i),
1359 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1361 write (iout,*) "SC restraints from templates"
1363 write(iout,'(i5,100(4f8.2,4x))') i,
1364 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1365 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1368 c -----------------------------------------------------------------
1371 c----------------------------------------------------------------------