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,1)
41 call readi(controlcard,'NSTART',nstart,0)
42 call readi(controlcard,'NEND',nend,0)
43 call reada(controlcard,'ECUT',ecut,10.0d0)
44 call reada(controlcard,'PROB',prob_limit,0.99d0)
45 write (iout,*) "Probability limit",prob_limit
46 lgrp=(index(controlcard,'LGRP').gt.0)
47 caonly=(index(controlcard,'CA_ONLY').gt.0)
48 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
49 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
50 call readi(controlcard,'IOPT',iopt,2)
51 lside = index(controlcard,"SIDE").gt.0
52 efree = index(controlcard,"EFREE").gt.0
53 call readi(controlcard,'NTEMP',nT,1)
54 write (iout,*) "nT",nT
55 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
56 write (iout,*) "nT",nT
57 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
59 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
61 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
62 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
63 lprint_int=index(controlcard,"PRINT_INT") .gt.0
64 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
65 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
66 write (iout,*) "with_dihed_constr ",with_dihed_constr,
67 & " CONSTR_DIST",constr_dist
69 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
70 write (iout,*) "with_homology_constr ",with_dihed_constr,
71 & " CONSTR_HOMOLOGY",constr_homology
77 c--------------------------------------------------------------------------
80 C Read molecular data.
84 include 'COMMON.IOUNITS'
87 include 'COMMON.INTERACT'
88 include 'COMMON.LOCAL'
89 include 'COMMON.NAMES'
90 include 'COMMON.CHAIN'
91 include 'COMMON.FFIELD'
92 include 'COMMON.SBRIDGE'
93 include 'COMMON.HEADER'
94 include 'COMMON.CONTROL'
95 include 'COMMON.CONTACTS'
96 include 'COMMON.TIME1'
97 include 'COMMON.TORCNSTR'
101 character*4 sequence(maxres)
102 character*800 weightcard
104 double precision x(maxvar)
105 integer itype_pdb(maxres)
108 write (iout,*) " MOLREAD: NRES",NRES
112 C Read weights of the subsequent energy terms.
113 call card_concat(weightcard)
114 call reada(weightcard,'WSC',wsc,1.0d0)
115 call reada(weightcard,'WLONG',wsc,wsc)
116 call reada(weightcard,'WSCP',wscp,1.0d0)
117 call reada(weightcard,'WELEC',welec,1.0D0)
118 call reada(weightcard,'WVDWPP',wvdwpp,welec)
119 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
120 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
121 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
122 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
123 call reada(weightcard,'WTURN3',wturn3,1.0D0)
124 call reada(weightcard,'WTURN4',wturn4,1.0D0)
125 call reada(weightcard,'WTURN6',wturn6,1.0D0)
126 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
127 call reada(weightcard,'WBOND',wbond,1.0D0)
128 call reada(weightcard,'WTOR',wtor,1.0D0)
129 call reada(weightcard,'WTORD',wtor_d,1.0D0)
130 call reada(weightcard,'WANG',wang,1.0D0)
131 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
132 call reada(weightcard,'SCAL14',scal14,0.4D0)
133 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
134 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
135 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
136 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
137 call reada(weightcard,"D0CM",d0cm,3.78d0)
138 call reada(weightcard,"AKCM",akcm,15.1d0)
139 call reada(weightcard,"AKTH",akth,11.0d0)
140 call reada(weightcard,"AKCT",akct,12.0d0)
141 call reada(weightcard,"V1SS",v1ss,-1.08d0)
142 call reada(weightcard,"V2SS",v2ss,7.61d0)
143 call reada(weightcard,"V3SS",v3ss,13.7d0)
144 call reada(weightcard,"EBR",ebr,-5.50D0)
146 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
147 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
148 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
149 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
150 if (index(weightcard,'SOFT').gt.0) ipot=6
151 C 12/1/95 Added weight for the multi-body term WCORR
152 call reada(weightcard,'WCORRH',wcorr,1.0D0)
155 dyn_ssbond_ij(i,j)=1.0d300
158 call reada(weightcard,"HT",Ht,0.0D0)
160 ss_depth=ebr/wsc-0.25*eps(1,1)
161 Ht=Ht/wsc-0.25*eps(1,1)
162 akcm=akcm*wstrain/wsc
163 akth=akth*wstrain/wsc
164 akct=akct*wstrain/wsc
165 v1ss=v1ss*wstrain/wsc
166 v2ss=v2ss*wstrain/wsc
167 v3ss=v3ss*wstrain/wsc
169 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
171 write (iout,'(/a)') "Disulfide bridge parameters:"
172 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
173 write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
174 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
175 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
176 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
178 write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
179 if (wcorr4.gt.0.0d0) wcorr=wcorr4
198 weights(22)=wdfa_dist
201 weights(25)=wdfa_beta
202 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
203 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
204 & wturn4,wturn6,wsccor,wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
205 10 format (/'Energy-term weights (unscaled):'//
206 & 'WSCC= ',f10.6,' (SC-SC)'/
207 & 'WSCP= ',f10.6,' (SC-p)'/
208 & 'WELEC= ',f10.6,' (p-p electr)'/
209 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
210 & 'WBOND= ',f10.6,' (stretching)'/
211 & 'WANG= ',f10.6,' (bending)'/
212 & 'WSCLOC= ',f10.6,' (SC local)'/
213 & 'WTOR= ',f10.6,' (torsional)'/
214 & 'WTORD= ',f10.6,' (double torsional)'/
215 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
216 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
217 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
218 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
219 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
220 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
221 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
222 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
223 & 'WSCCOR= ',f10.6,' (SC-backbone torsional correalations)'/
224 & 'WDFAD= ',f10.6,' (DFA distance)'/
225 & 'WDFAT= ',f10.6,' (DFA torsional)'/
226 & 'WDFAN= ',f10.6,' (DFA neighbors)'/
227 & 'WDFAB= ',f10.6,' (DFA beta)'/)
228 if (wcorr4.gt.0.0d0) then
229 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
230 & 'between contact pairs of peptide groups'
231 write (iout,'(2(a,f5.3/))')
232 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
233 & 'Range of quenching the correlation terms:',2*delt_corr
234 else if (wcorr.gt.0.0d0) then
235 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
236 & 'between contact pairs of peptide groups'
238 write (iout,'(a,f8.3)')
239 & 'Scaling factor of 1,4 SC-p interactions:',scal14
240 write (iout,'(a,f8.3)')
241 & 'General scaling factor of SC-p interactions:',scalscp
242 r0_corr=cutoff_corr-delt_corr
244 aad(i,1)=scalscp*aad(i,1)
245 aad(i,2)=scalscp*aad(i,2)
246 bad(i,1)=scalscp*bad(i,1)
247 bad(i,2)=scalscp*bad(i,2)
251 print *,'indpdb=',indpdb,' pdbref=',pdbref
253 C Read sequence if not taken from the pdb file.
254 if (iscode.gt.0) then
255 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
257 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
259 C Convert sequence to numeric code
261 itype(i)=rescode(i,sequence(i),iscode)
264 print '(20i4)',(itype(i),i=1,nres)
268 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
270 if (itype(i).eq.21) then
274 else if (itype(i+1).ne.20) then
276 else if (itype(i).ne.20) then
283 write (iout,*) "ITEL"
285 write (iout,*) i,itype(i),itel(i)
288 print *,'Call Read_Bridge.'
290 if (with_dihed_constr) then
292 read (inp,*) ndih_constr
293 if (ndih_constr.gt.0) then
295 write (iout,*) 'FTORS',ftors
296 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
298 & 'There are',ndih_constr,' constraints on phi angles.'
300 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
303 phi0(i)=deg2rad*phi0(i)
304 drange(i)=deg2rad*drange(i)
312 print *,'NNT=',NNT,' NCT=',NCT
313 if (itype(1).eq.21) nnt=2
314 if (itype(nres).eq.21) nct=nct-1
315 if (nstart.lt.nnt) nstart=nnt
316 if (nend.gt.nct .or. nend.eq.0) nend=nct
317 write (iout,*) "nstart",nstart," nend",nend
320 C Juyong:READ init_vars
321 C Initialize variables!
322 C Juyong:READ read_info
323 C READ fragment information!!
324 C both routines should be in dfa.F file!!
326 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
327 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
328 write (iout,*) "Calling init_dfa_vars"
331 write (iout,*) 'init_dfa_vars finished!'
334 write (iout,*) 'read_dfa_info finished!'
338 if (constr_homology.gt.0) then
339 call read_constr_homology
343 c read(inp,'(a)') pdbfile
344 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
345 c open(ipdbin,file=pdbfile,status='old',err=33)
347 c 33 write (iout,'(a)') 'Error opening PDB file.'
350 c print *,'Begin reading pdb data'
352 c print *,'Finished reading pdb data'
353 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
355 c itype_pdb(i)=itype(i)
358 c write (iout,'(a,i3)') 'nsup=',nsup
360 c if (nsup.le.(nct-nnt+1)) then
361 c do i=0,nct-nnt+1-nsup
362 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
368 c & 'Error - sequences to be superposed do not match.'
371 c do i=0,nsup-(nct-nnt+1)
372 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
374 c nstart_sup=nstart_sup+i
380 c & 'Error - sequences to be superposed do not match.'
383 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
384 c & ' nstart_seq=',nstart_seq
388 write (iout,*) "molread: REFSTR",refstr
390 if (.not.pdbref) then
391 call read_angles(inp,*38)
393 38 write (iout,'(a)') 'Error reading reference structure.'
395 call mp_stopall(Error_Msg)
397 stop 'Error reading reference structure'
409 call contact(.true.,ncont_ref,icont_ref)
411 c Read distance restraints
412 if (constr_dist.gt.0) then
413 call read_dist_constr
418 c-----------------------------------------------------------------------------
419 logical function seq_comp(itypea,itypeb,length)
421 integer length,itypea(length),itypeb(length)
424 if (itypea(i).ne.itypeb(i)) then
432 c-----------------------------------------------------------------------------
433 subroutine read_bridge
434 C Read information about disulfide bridges.
437 include 'COMMON.IOUNITS'
440 include 'COMMON.INTERACT'
441 include 'COMMON.LOCAL'
442 include 'COMMON.NAMES'
443 include 'COMMON.CHAIN'
444 include 'COMMON.FFIELD'
445 include 'COMMON.SBRIDGE'
446 include 'COMMON.HEADER'
447 include 'COMMON.CONTROL'
448 include 'COMMON.TIME1'
450 include 'COMMON.INFO'
453 C Read bridging residues.
454 read (inp,*) ns,(iss(i),i=1,ns)
455 write(iout,*)'ns=',ns
456 C Check whether the specified bridging residues are cystines.
458 if (itype(iss(i)).ne.1) then
459 write (iout,'(2a,i3,a)')
460 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
461 & ' can form a disulfide bridge?!!!'
462 write (*,'(2a,i3,a)')
463 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
464 & ' can form a disulfide bridge?!!!'
466 call mp_stopall(error_msg)
472 C Read preformed bridges.
474 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
477 C Check if the residues involved in bridges are in the specified list of
481 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
482 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
483 write (iout,'(a,i3,a)') 'Disulfide pair',i,
484 & ' contains residues present in other pairs.'
485 write (*,'(a,i3,a)') 'Disulfide pair',i,
486 & ' contains residues present in other pairs.'
488 call mp_stopall(error_msg)
495 if (ihpb(i).eq.iss(j)) goto 10
497 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
500 if (jhpb(i).eq.iss(j)) goto 20
502 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
513 if (ns.gt.0.and.dyn_ss) then
517 forcon(i-nss)=forcon(i)
524 dyn_ss_mask(iss(i))=.true.
525 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
528 print *, "Leaving brigde read"
531 c----------------------------------------------------------------------------
532 subroutine read_angles(kanal,*)
537 include 'COMMON.CHAIN'
538 include 'COMMON.IOUNITS'
540 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
541 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
542 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
543 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
545 theta(i)=deg2rad*theta(i)
546 phi(i)=deg2rad*phi(i)
547 alph(i)=deg2rad*alph(i)
548 omeg(i)=deg2rad*omeg(i)
553 c----------------------------------------------------------------------------
554 subroutine reada(rekord,lancuch,wartosc,default)
556 character*(*) rekord,lancuch
557 double precision wartosc,default
560 iread=index(rekord,lancuch)
565 iread=iread+ilen(lancuch)+1
566 read (rekord(iread:),*) wartosc
569 c----------------------------------------------------------------------------
570 subroutine multreada(rekord,lancuch,tablica,dim,default)
573 double precision tablica(dim),default
574 character*(*) rekord,lancuch
580 iread=index(rekord,lancuch)
581 if (iread.eq.0) return
582 iread=iread+ilen(lancuch)+1
583 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
586 c----------------------------------------------------------------------------
587 subroutine readi(rekord,lancuch,wartosc,default)
589 character*(*) rekord,lancuch
590 integer wartosc,default
593 iread=index(rekord,lancuch)
598 iread=iread+ilen(lancuch)+1
599 read (rekord(iread:),*) wartosc
602 c----------------------------------------------------------------------------
603 subroutine multreadi(rekord,lancuch,tablica,dim,default)
606 integer tablica(dim),default
607 character*(*) rekord,lancuch
614 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
615 if (iread.eq.0) return
616 iread=iread+ilen(lancuch)+1
617 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
620 c----------------------------------------------------------------------------
621 subroutine card_concat(card)
623 include 'COMMON.IOUNITS'
625 character*80 karta,ucase
627 read (inp,'(a)') karta
630 do while (karta(80:80).eq.'&')
631 card=card(:ilen(card)+1)//karta(:79)
632 read (inp,'(a)') karta
635 card=card(:ilen(card)+1)//karta
638 c----------------------------------------------------------------------------
647 include 'COMMON.IOUNITS'
648 include 'COMMON.CONTROL'
649 integer lenpre,lenpot,ilen
651 character*16 cformat,cprint
653 integer lenint,lenout
654 call getenv('INPUT',prefix)
655 call getenv('OUTPUT',prefout)
656 call getenv('INTIN',prefintin)
657 call getenv('COORD',cformat)
658 call getenv('PRINTCOOR',cprint)
659 call getenv('SCRATCHDIR',scratchdir)
662 if (index(ucase(cformat),'CX').gt.0) then
669 lenint=ilen(prefintin)
670 C Get the names and open the input files
671 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
673 write (liczba,'(bz,i3.3)') me
674 outname=prefout(:lenout)//'_clust.out_'//liczba
676 outname=prefout(:lenout)//'_clust.out'
679 intinname=prefintin(:lenint)//'.bx'
680 else if (from_cx) then
681 intinname=prefintin(:lenint)//'.cx'
683 intinname=prefintin(:lenint)//'.int'
685 rmsname=prefintin(:lenint)//'.rms'
686 open (jplot,file=prefout(:ilen(prefout))//'.tex',
688 open (jrms,file=rmsname,status='unknown')
689 open(iout,file=outname,status='unknown')
690 C Get parameter filenames and open the parameter files.
691 call getenv('BONDPAR',bondname)
692 open (ibond,file=bondname,status='old')
693 call getenv('THETPAR',thetname)
694 open (ithep,file=thetname,status='old')
695 call getenv('ROTPAR',rotname)
696 open (irotam,file=rotname,status='old')
697 call getenv('TORPAR',torname)
698 open (itorp,file=torname,status='old')
699 call getenv('TORDPAR',tordname)
700 open (itordp,file=tordname,status='old')
701 call getenv('FOURIER',fouriername)
702 open (ifourier,file=fouriername,status='old')
703 call getenv('ELEPAR',elename)
704 open (ielep,file=elename,status='old')
705 call getenv('SIDEPAR',sidename)
706 open (isidep,file=sidename,status='old')
707 call getenv('SIDEP',sidepname)
708 open (isidep1,file=sidepname,status="old")
709 call getenv('SCCORPAR',sccorname)
710 open (isccor,file=sccorname,status="old")
713 C 8/9/01 In the newest version SCp interaction constants are read from a file
714 C Use -DOLDSCP to use hard-coded constants instead.
716 call getenv('SCPPAR',scpname)
717 open (iscpp,file=scpname,status='old')
721 c-------------------------------------------------------------------------------
722 subroutine read_dist_constr
723 implicit real*8 (a-h,o-z)
725 include 'COMMON.CONTROL'
726 include 'COMMON.CHAIN'
727 include 'COMMON.IOUNITS'
728 include 'COMMON.SBRIDGE'
729 integer ifrag_(2,100),ipair_(2,100)
730 double precision wfrag_(100),wpair_(100)
731 character*500 controlcard
732 c write (iout,*) "Calling read_dist_constr"
733 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
734 call card_concat(controlcard)
736 c write (iout,'(a)') controlcard
737 call readi(controlcard,"NFRAG",nfrag_,0)
738 call readi(controlcard,"NPAIR",npair_,0)
739 call readi(controlcard,"NDIST",ndist_,0)
740 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
741 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
742 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
743 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
744 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
745 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
746 write (iout,*) "IFRAG"
748 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
750 write (iout,*) "IPAIR"
752 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
755 if (.not.refstr .and. nfrag_.gt.0) then
757 & "ERROR: no reference structure to compute distance restraints"
759 & "Restraints must be specified explicitly (NDIST=number)"
762 if (nfrag_.lt.2 .and. npair_.gt.0) then
763 write (iout,*) "ERROR: Less than 2 fragments specified",
764 & " but distance restraints between pairs requested"
769 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
770 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
771 & ifrag_(2,i)=nstart_sup+nsup-1
772 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
774 if (wfrag_(i).gt.0.0d0) then
775 do j=ifrag_(1,i),ifrag_(2,i)-1
777 write (iout,*) "j",j," k",k
779 if (constr_dist.eq.1) then
784 forcon(nhpb)=wfrag_(i)
785 else if (constr_dist.eq.2) then
786 if (ddjk.le.dist_cut) then
791 forcon(nhpb)=wfrag_(i)
798 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
800 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
801 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
807 if (wpair_(i).gt.0.0d0) then
815 do j=ifrag_(1,ii),ifrag_(2,ii)
816 do k=ifrag_(1,jj),ifrag_(2,jj)
820 forcon(nhpb)=wpair_(i)
822 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
823 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
829 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
830 & ibecarb(i),forcon(nhpb+1)
831 if (forcon(nhpb+1).gt.0.0d0) then
833 if (ibecarb(i).gt.0) then
837 if (dhpb(nhpb).eq.0.0d0)
838 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
842 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
843 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
849 c====-------------------------------------------------------------------
850 subroutine read_constr_homology
856 include 'COMMON.SETUP'
857 include 'COMMON.CONTROL'
858 include 'COMMON.CHAIN'
859 include 'COMMON.IOUNITS'
861 include 'COMMON.INTERACT'
862 include 'COMMON.HOMRESTR'
867 c include 'include_unres/COMMON.VAR'
870 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
872 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
873 c & sigma_odl_temp(maxres,maxres,max_template)
875 character*24 model_ki_dist, model_ki_angle
876 character*500 controlcard
877 integer ki, i, j, k, l
878 logical lprn /.true./
880 c FP - Nov. 2014 Temporary specifications for new vars
882 double precision rescore_tmp,x12,y12,z12
883 double precision, dimension (max_template,maxres) :: rescore
884 character*24 tpl_k_rescore
885 c -----------------------------------------------------------------
886 c Reading multiple PDB ref structures and calculation of retraints
887 c not using pre-computed ones stored in files model_ki_{dist,angle}
889 c -----------------------------------------------------------------
892 c Alternative: reading from input
893 write (iout,*) "BEGIN READ HOMOLOGY INFO"
895 call card_concat(controlcard)
896 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
897 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
898 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
899 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
900 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
902 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
903 if (homol_nset.gt.1)then
904 call readi(controlcard,"ISET",iset,homol_nset)
905 call card_concat(controlcard)
906 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
912 write(iout,*) "read_constr_homology iset",iset
913 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
917 cd write (iout,*) "nnt",nnt," nct",nct
929 c Reading HM global scores (prob not required)
931 c open (4,file="HMscore")
932 c do k=1,constr_homology
933 c read (4,*,end=521) hmscore_tmp
934 c hmscore(k)=hmscore_tmp ! Another transformation can be used
935 c write(*,*) "Model", k, ":", hmscore(k)
939 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
941 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
942 do k=1,constr_homology
944 read(inp,'(a)') pdbfile
945 write (iout,*) "k ",k," pdbfile ",pdbfile
946 c Next stament causes error upon compilation (?)
947 c if(me.eq.king.or. .not. out1file)
948 c write (iout,'(2a)') 'PDB data will be read from file ',
949 c & pdbfile(:ilen(pdbfile))
950 open(ipdbin,file=pdbfile,status='old',err=33)
952 33 write (iout,'(a)') 'Error opening PDB file.'
955 c print *,'Begin reading pdb data'
957 c Files containing res sim or local scores (former containing sigmas)
960 write(kic2,'(bz,i2.2)') k
962 tpl_k_rescore="template"//kic2//".sco"
963 c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
964 c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
965 c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
966 c tpl_k_sigma_d="template"//kic2//".sigma_d"
977 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
978 & (crefjlee(j,i+nres),j=1,3)
981 write (iout,*) "READ HOMOLOGY INFO"
982 write (iout,*) "read_constr_homology x: after reading pdb file"
983 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
984 write (iout,*) "waga_dist",waga_dist
985 write (iout,*) "waga_angle",waga_angle
986 write (iout,*) "waga_theta",waga_theta
987 write (iout,*) "waga_d",waga_d
988 write (iout,*) "dist_cut",dist_cut
992 c Distance restraints
995 C Copy the coordinates from reference coordinates (?)
999 c write (iout,*) "c(",j,i,") =",c(j,i)
1003 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1005 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1006 open (ientin,file=tpl_k_rescore,status='old')
1007 do irec=1,maxdim ! loop for reading res sim
1009 rescore(k,irec)=0.0d0
1012 read (ientin,*,end=1401) rescore_tmp
1013 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1014 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1015 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1020 c open (ientin,file=tpl_k_sigma_odl,status='old')
1021 c do irec=1,maxdim ! loop for reading sigma_odl
1022 c read (ientin,*,end=1401) i, j,
1023 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
1024 c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
1025 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
1029 if (waga_dist.ne.0.0d0) then
1031 do i = nnt,nct-2 ! right? without parallel.
1032 do j=i+2,nct ! right?
1033 c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
1034 c do j=i+2,nres ! ibid
1035 c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
1036 c do j=i+2,nct ! ibid
1038 c write (iout,*) "k",k
1039 c write (iout,*) "i",i," j",j," constr_homology",
1044 c Attempt to replace dist(i,j) by its definition in ...
1049 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1052 c odl(k,ii)=dist(i,j)
1053 c write (iout,*) "dist(",i,j,") =",dist(i,j)
1054 c write (iout,*) "distal = ",distal
1055 c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
1056 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1057 c & "rescore(",k,j,") =",rescore(k,j)
1059 c Calculation of sigma from res sim
1061 c if (odl(k,ii).le.6.0d0) then
1062 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
1063 c Other functional forms possible depending on odl(k,ii), eg.
1065 if (odl(k,ii).le.dist_cut) then
1066 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
1067 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
1070 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
1071 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1073 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
1074 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1076 c Following expr replaced by a positive exp argument
1077 c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1078 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
1080 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
1081 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
1084 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
1085 c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
1087 c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
1088 c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
1090 c read (ientin,*) sigma_odl(k,ii) ! 1st variant
1093 c if (constr_homology.gt.0) call homology_partition
1096 c Theta, dihedral and SC retraints
1098 if (waga_angle.gt.0.0d0) then
1099 c open (ientin,file=tpl_k_sigma_dih,status='old')
1100 c do irec=1,maxres-3 ! loop for reading sigma_dih
1101 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1102 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1103 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1104 c & sigma_dih(k,i+nnt-1)
1108 do i = nnt+3,nct ! right? without parallel.
1109 c do i=1,nres ! alternative for bounds acc to readpdb?
1110 c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
1111 c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
1112 dih(k,i)=phiref(i) ! right?
1113 c read (ientin,*) sigma_dih(k,i) ! original variant
1114 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1115 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1116 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1117 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1118 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1120 sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
1121 & rescore(k,i-2)+rescore(k,i-3) ! right expression ?
1123 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1124 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1125 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1126 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1127 c Instead of res sim other local measure of b/b str reliability possible
1128 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1129 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1130 if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
1131 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
1135 if (waga_theta.gt.0.0d0) then
1136 c open (ientin,file=tpl_k_sigma_theta,status='old')
1137 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1138 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1139 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1140 c & sigma_theta(k,i+nnt-1)
1145 do i = nnt+2,nct ! right? without parallel.
1146 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1147 c do i=ithet_start,ithet_end ! with FG parallel.
1148 thetatpl(k,i)=thetaref(i)
1149 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1150 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1151 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1152 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1153 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1154 sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
1155 & rescore(k,i-2) ! right expression ?
1156 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1158 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1159 c rescore(k,i-2) ! right expression ?
1160 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1161 if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
1165 if (waga_d.gt.0.0d0) then
1166 c open (ientin,file=tpl_k_sigma_d,status='old')
1167 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1168 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1169 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1170 c & sigma_d(k,i+nnt-1)
1175 do i = nnt,nct ! right? without parallel.
1176 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1177 c do i=loc_start,loc_end ! with FG parallel.
1178 if (itype(i).eq.10) goto 1 ! right?
1182 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1183 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1184 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1185 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1186 sigma_d(k,i)=rescore(k,i) ! right expression ?
1187 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1189 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1190 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1191 c read (ientin,*) sigma_d(k,i) ! 1st variant
1192 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1198 if (waga_dist.ne.0.0d0) lim_odl=ii
1199 if (constr_homology.gt.0) call homology_partition
1200 if (constr_homology.gt.0) call init_int_table
1201 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1202 cd & "lim_xx=",lim_xx
1203 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1204 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1208 if (.not.lprn) return
1209 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1210 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1211 write (iout,*) "Distance restraints from templates"
1213 write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
1214 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
1216 write (iout,*) "Dihedral angle restraints from templates"
1218 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
1219 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1221 write (iout,*) "Virtual-bond angle restraints from templates"
1222 do i=nnt+2,lim_theta
1223 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
1224 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1226 write (iout,*) "SC restraints from templates"
1228 write(iout,'(i5,10(4f8.2,4x))') i,
1229 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1230 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1233 c -----------------------------------------------------------------