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 include 'COMMON.INTERACT'
18 include "COMMON.SPLITELE"
19 include 'COMMON.SHIELD'
20 character*320 controlcard,ucase
24 integer i,i1,i2,it1,it2
26 read (INP,'(a80)') titel
27 call card_concat(controlcard)
29 call readi(controlcard,'NRES',nres,0)
30 call readi(controlcard,'RESCALE',rescale_mode,2)
31 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
32 write (iout,*) "DISTCHAINMAX",distchainmax
33 C Reading the dimensions of box in x,y,z coordinates
34 call reada(controlcard,'BOXX',boxxsize,100.0d0)
35 call reada(controlcard,'BOXY',boxysize,100.0d0)
36 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
37 c Cutoff range for interactions
38 call reada(controlcard,"R_CUT",r_cut,15.0d0)
39 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
40 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
41 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
42 if (lipthick.gt.0.0d0) then
43 bordliptop=(boxzsize+lipthick)/2.0
44 bordlipbot=bordliptop-lipthick
46 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
47 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
48 buflipbot=bordlipbot+lipbufthick
49 bufliptop=bordliptop-lipbufthick
50 if ((lipbufthick*2.0d0).gt.lipthick)
51 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
53 write(iout,*) "bordliptop=",bordliptop
54 write(iout,*) "bordlipbot=",bordlipbot
55 write(iout,*) "bufliptop=",bufliptop
56 write(iout,*) "buflipbot=",buflipbot
58 call readi(controlcard,'SHIELD',shield_mode,0)
59 write (iout,*) "SHIELD MODE",shield_mode
60 if (shield_mode.gt.0) then
62 C VSolvSphere the volume of solving sphere
64 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
65 C there will be no distinction between proline peptide group and normal peptide
66 C group in case of shielding parameters
67 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
68 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
69 write (iout,*) VSolvSphere,VSolvSphere_div
70 C long axis of side chain
72 long_r_sidechain(i)=vbldsc0(1,i)
73 short_r_sidechain(i)=sigma0(i)
77 call readi(controlcard,'PDBOUT',outpdb,0)
78 call readi(controlcard,'MOL2OUT',outmol2,0)
79 refstr=(index(controlcard,'REFSTR').gt.0)
80 write (iout,*) "REFSTR",refstr
81 pdbref=(index(controlcard,'PDBREF').gt.0)
82 iscode=index(controlcard,'ONE_LETTER')
83 tree=(index(controlcard,'MAKE_TREE').gt.0)
84 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
85 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
86 write (iout,*) "with_dihed_constr ",with_dihed_constr,
87 & " CONSTR_DIST",constr_dist
88 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
89 write (iout,*) "with_theta_constr ",with_theta_constr
91 min_var=(index(controlcard,'MINVAR').gt.0)
92 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
93 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
94 call readi(controlcard,'NCUT',ncut,0)
95 call readi(controlcard,'NCLUST',nclust,5)
96 call readi(controlcard,'SYM',symetr,1)
97 write (iout,*) 'sym', symetr
98 call readi(controlcard,'NSTART',nstart,0)
99 call readi(controlcard,'NEND',nend,0)
100 call reada(controlcard,'ECUT',ecut,10.0d0)
101 call reada(controlcard,'PROB',prob_limit,0.99d0)
102 write (iout,*) "Probability limit",prob_limit
103 lgrp=(index(controlcard,'LGRP').gt.0)
104 caonly=(index(controlcard,'CA_ONLY').gt.0)
105 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
107 & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
108 call readi(controlcard,'IOPT',iopt,2)
109 lside = index(controlcard,"SIDE").gt.0
110 efree = index(controlcard,"EFREE").gt.0
111 call readi(controlcard,'NTEMP',nT,1)
112 write (iout,*) "nT",nT
113 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
114 write (iout,*) "nT",nT
115 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
117 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
119 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
120 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
121 lprint_int=index(controlcard,"PRINT_INT") .gt.0
122 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
123 write (iout,*) "with_homology_constr ",with_dihed_constr,
124 & " CONSTR_HOMOLOGY",constr_homology
125 print_homology_restraints=
126 & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
127 print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
128 print_homology_models=
129 & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
133 c--------------------------------------------------------------------------
136 C Read molecular data.
140 include 'COMMON.IOUNITS'
143 include 'COMMON.INTERACT'
144 include 'COMMON.LOCAL'
145 include 'COMMON.NAMES'
146 include 'COMMON.CHAIN'
147 include 'COMMON.FFIELD'
148 include 'COMMON.SBRIDGE'
149 include 'COMMON.HEADER'
150 include 'COMMON.CONTROL'
151 include 'COMMON.CONTACTS'
152 include 'COMMON.TIME1'
153 include 'COMMON.TORCNSTR'
154 include 'COMMON.SHIELD'
156 include 'COMMON.INFO'
158 character*4 sequence(maxres)
159 character*800 weightcard
161 double precision x(maxvar)
162 integer itype_pdb(maxres)
164 integer i,j,kkk,i1,i2,it1,it2
168 C Read weights of the subsequent energy terms.
169 call card_concat(weightcard)
170 write(iout,*) weightcard
171 C call reada(weightcard,'WSC',wsc,1.0d0)
173 call reada(weightcard,'WLONG',wsc,wsc)
174 call reada(weightcard,'WSCP',wscp,1.0d0)
175 call reada(weightcard,'WELEC',welec,1.0D0)
176 call reada(weightcard,'WVDWPP',wvdwpp,welec)
177 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
178 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
179 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
180 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
181 call reada(weightcard,'WTURN3',wturn3,1.0D0)
182 call reada(weightcard,'WTURN4',wturn4,1.0D0)
183 call reada(weightcard,'WTURN6',wturn6,1.0D0)
184 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
185 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
186 call reada(weightcard,'WBOND',wbond,1.0D0)
187 call reada(weightcard,'WTOR',wtor,1.0D0)
188 call reada(weightcard,'WTORD',wtor_d,1.0D0)
189 call reada(weightcard,'WANG',wang,1.0D0)
190 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
191 call reada(weightcard,'SCAL14',scal14,0.4D0)
192 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
193 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
194 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
195 if (index(weightcard,'SOFT').gt.0) ipot=6
196 call reada(weightcard,"D0CM",d0cm,3.78d0)
197 call reada(weightcard,"AKCM",akcm,15.1d0)
198 call reada(weightcard,"AKTH",akth,11.0d0)
199 call reada(weightcard,"AKCT",akct,12.0d0)
200 call reada(weightcard,"V1SS",v1ss,-1.08d0)
201 call reada(weightcard,"V2SS",v2ss,7.61d0)
202 call reada(weightcard,"V3SS",v3ss,13.7d0)
203 call reada(weightcard,"EBR",ebr,-5.50D0)
204 call reada(weightcard,'WSHIELD',wshield,1.0d0)
205 write(iout,*) 'WSHIELD',wshield
206 call reada(weightcard,'WLT',wliptran,0.0D0)
207 call reada(weightcard,"ATRISS",atriss,0.301D0)
208 call reada(weightcard,"BTRISS",btriss,0.021D0)
209 call reada(weightcard,"CTRISS",ctriss,1.001D0)
210 call reada(weightcard,"DTRISS",dtriss,1.001D0)
211 write (iout,*) "ATRISS=", atriss
212 write (iout,*) "BTRISS=", btriss
213 write (iout,*) "CTRISS=", ctriss
214 write (iout,*) "DTRISS=", dtriss
215 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
217 dyn_ss_mask(i)=.false.
221 dyn_ssbond_ij(i,j)=1.0d300
224 call reada(weightcard,"HT",Ht,0.0D0)
226 ss_depth=ebr/wsc-0.25*eps(1,1)
227 Ht=Ht/wsc-0.25*eps(1,1)
228 akcm=akcm*wstrain/wsc
229 akth=akth*wstrain/wsc
230 akct=akct*wstrain/wsc
231 v1ss=v1ss*wstrain/wsc
232 v2ss=v2ss*wstrain/wsc
233 v3ss=v3ss*wstrain/wsc
235 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
237 write (iout,'(/a)') "Disulfide bridge parameters:"
238 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
239 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
240 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
241 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
244 C 12/1/95 Added weight for the multi-body term WCORR
245 call reada(weightcard,'WCORRH',wcorr,1.0D0)
246 if (wcorr4.gt.0.0d0) wcorr=wcorr4
266 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
267 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
268 & wturn4,wturn6,wsccor
269 10 format (/'Energy-term weights (unscaled):'//
270 & 'WSCC= ',f10.6,' (SC-SC)'/
271 & 'WSCP= ',f10.6,' (SC-p)'/
272 & 'WELEC= ',f10.6,' (p-p electr)'/
273 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
274 & 'WBOND= ',f10.6,' (stretching)'/
275 & 'WANG= ',f10.6,' (bending)'/
276 & 'WSCLOC= ',f10.6,' (SC local)'/
277 & 'WTOR= ',f10.6,' (torsional)'/
278 & 'WTORD= ',f10.6,' (double torsional)'/
279 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
280 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
281 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
282 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
283 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
284 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
285 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
286 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
287 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
289 if (wcorr4.gt.0.0d0) then
290 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
291 & 'between contact pairs of peptide groups'
292 write (iout,'(2(a,f5.3/))')
293 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
294 & 'Range of quenching the correlation terms:',2*delt_corr
295 else if (wcorr.gt.0.0d0) then
296 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
297 & 'between contact pairs of peptide groups'
299 write (iout,'(a,f8.3)')
300 & 'Scaling factor of 1,4 SC-p interactions:',scal14
301 write (iout,'(a,f8.3)')
302 & 'General scaling factor of SC-p interactions:',scalscp
303 r0_corr=cutoff_corr-delt_corr
305 aad(i,1)=scalscp*aad(i,1)
306 aad(i,2)=scalscp*aad(i,2)
307 bad(i,1)=scalscp*bad(i,1)
308 bad(i,2)=scalscp*bad(i,2)
312 print *,'indpdb=',indpdb,' pdbref=',pdbref
314 C Read sequence if not taken from the pdb file.
315 if (iscode.gt.0) then
316 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
318 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
320 C Convert sequence to numeric code
322 itype(i)=rescode(i,sequence(i),iscode)
325 print '(20i4)',(itype(i),i=1,nres)
329 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
331 if (itype(i).eq.ntyp1) then
335 else if (iabs(itype(i+1)).ne.20) then
337 else if (iabs(itype(i)).ne.20) then
344 write (iout,*) "ITEL"
346 write (iout,*) i,itype(i),itel(i)
349 print *,'Call Read_Bridge.'
351 C this fragment reads diheadral constrains
352 if (with_dihed_constr) then
354 read (inp,*) ndih_constr
355 if (ndih_constr.gt.0) then
357 C write (iout,*) 'FTORS',ftors
358 C ftors is the force constant for torsional quartic constrains
359 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
362 & 'There are',ndih_constr,' constraints on phi angles.'
364 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
368 phi0(i)=deg2rad*phi0(i)
369 drange(i)=deg2rad*drange(i)
371 endif ! endif ndif_constr.gt.0
372 endif ! with_dihed_constr
373 if (with_theta_constr) then
374 C with_theta_constr is keyword allowing for occurance of theta constrains
375 read (inp,*) ntheta_constr
376 C ntheta_constr is the number of theta constrains
377 if (ntheta_constr.gt.0) then
379 read (inp,*) (itheta_constr(i),theta_constr0(i),
380 & theta_drange(i),for_thet_constr(i),
382 C the above code reads from 1 to ntheta_constr
383 C itheta_constr(i) residue i for which is theta_constr
384 C theta_constr0 the global minimum value
385 C theta_drange is range for which there is no energy penalty
386 C for_thet_constr is the force constant for quartic energy penalty
388 C if(me.eq.king.or..not.out1file)then
390 & 'There are',ntheta_constr,' constraints on phi angles.'
392 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
398 theta_constr0(i)=deg2rad*theta_constr0(i)
399 theta_drange(i)=deg2rad*theta_drange(i)
401 C if(me.eq.king.or..not.out1file)
402 C & write (iout,*) 'FTORS',ftors
403 C do i=1,ntheta_constr
404 C ii = itheta_constr(i)
405 C thetabound(1,ii) = phi0(i)-drange(i)
406 C thetabound(2,ii) = phi0(i)+drange(i)
408 endif ! ntheta_constr.gt.0
409 endif! with_theta_constr
413 print *,'NNT=',NNT,' NCT=',NCT
414 if (itype(1).eq.ntyp1) nnt=2
415 if (itype(nres).eq.ntyp1) nct=nct-1
416 if (nstart.lt.nnt) nstart=nnt
417 if (nend.gt.nct .or. nend.eq.0) nend=nct
418 write (iout,*) "nstart",nstart," nend",nend
420 if (constr_homology.gt.0) then
421 call read_constr_homology(print_homology_restraints)
425 c read(inp,'(a)') pdbfile
426 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
427 c open(ipdbin,file=pdbfile,status='old',err=33)
429 c 33 write (iout,'(a)') 'Error opening PDB file.'
432 c print *,'Begin reading pdb data'
434 c print *,'Finished reading pdb data'
435 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
437 c itype_pdb(i)=itype(i)
440 c write (iout,'(a,i3)') 'nsup=',nsup
442 c if (nsup.le.(nct-nnt+1)) then
443 c do i=0,nct-nnt+1-nsup
444 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
450 c & 'Error - sequences to be superposed do not match.'
453 c do i=0,nsup-(nct-nnt+1)
454 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
456 c nstart_sup=nstart_sup+i
462 c & 'Error - sequences to be superposed do not match.'
465 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
466 c & ' nstart_seq=',nstart_seq
470 write (iout,*) "molread: REFSTR",refstr
472 if (.not.pdbref) then
473 call read_angles(inp,*38)
475 38 write (iout,'(a)') 'Error reading reference structure.'
477 call mp_stopall(Error_Msg)
479 stop 'Error reading reference structure'
492 call contact(.true.,ncont_ref,icont_ref)
495 C write (iout,'(/a,i3,a)')
496 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
497 write (iout,'(20i4)') (iss(i),i=1,ns)
499 write(iout,*)"Running with dynamic disulfide-bond formation"
501 write (iout,'(/a/)') 'Pre-formed links are:'
507 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
508 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
514 if (ns.gt.0.and.dyn_ss) then
518 forcon(i-nss)=forcon(i)
525 dyn_ss_mask(iss(i))=.true.
528 c Read distance restraints
529 if (constr_dist.gt.0) then
530 call read_dist_constr
535 c-----------------------------------------------------------------------------
536 logical function seq_comp(itypea,itypeb,length)
538 integer length,itypea(length),itypeb(length)
541 if (itypea(i).ne.itypeb(i)) then
549 c-----------------------------------------------------------------------------
550 subroutine read_bridge
551 C Read information about disulfide bridges.
554 include 'COMMON.IOUNITS'
557 include 'COMMON.INTERACT'
558 include 'COMMON.LOCAL'
559 include 'COMMON.NAMES'
560 include 'COMMON.CHAIN'
561 include 'COMMON.FFIELD'
562 include 'COMMON.SBRIDGE'
563 include 'COMMON.HEADER'
564 include 'COMMON.CONTROL'
565 include 'COMMON.TIME1'
567 include 'COMMON.INFO'
570 C Read bridging residues.
571 read (inp,*) ns,(iss(i),i=1,ns)
573 C Check whether the specified bridging residues are cystines.
575 if (itype(iss(i)).ne.1) then
576 write (iout,'(2a,i3,a)')
577 & 'Do you REALLY think that the residue ',
578 & restyp(itype(iss(i))),i,
579 & ' can form a disulfide bridge?!!!'
580 write (*,'(2a,i3,a)')
581 & 'Do you REALLY think that the residue ',
582 & restyp(itype(iss(i))),i,
583 & ' can form a disulfide bridge?!!!'
585 call mp_stopall(error_msg)
591 C Read preformed bridges.
593 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
596 C Check if the residues involved in bridges are in the specified list of
600 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
601 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
602 write (iout,'(a,i3,a)') 'Disulfide pair',i,
603 & ' contains residues present in other pairs.'
604 write (*,'(a,i3,a)') 'Disulfide pair',i,
605 & ' contains residues present in other pairs.'
607 call mp_stopall(error_msg)
614 if (ihpb(i).eq.iss(j)) goto 10
616 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
619 if (jhpb(i).eq.iss(j)) goto 20
621 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
634 c----------------------------------------------------------------------------
635 subroutine read_angles(kanal,*)
640 include 'COMMON.CHAIN'
641 include 'COMMON.IOUNITS'
643 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
644 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
645 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
646 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
648 theta(i)=deg2rad*theta(i)
649 phi(i)=deg2rad*phi(i)
650 alph(i)=deg2rad*alph(i)
651 omeg(i)=deg2rad*omeg(i)
656 c----------------------------------------------------------------------------
657 subroutine reada(rekord,lancuch,wartosc,default)
659 character*(*) rekord,lancuch
660 double precision wartosc,default
663 iread=index(rekord,lancuch)
668 iread=iread+ilen(lancuch)+1
669 read (rekord(iread:),*) wartosc
672 c----------------------------------------------------------------------------
673 subroutine multreada(rekord,lancuch,tablica,dim,default)
676 double precision tablica(dim),default
677 character*(*) rekord,lancuch
683 iread=index(rekord,lancuch)
684 if (iread.eq.0) return
685 iread=iread+ilen(lancuch)+1
686 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
689 c----------------------------------------------------------------------------
690 subroutine readi(rekord,lancuch,wartosc,default)
692 character*(*) rekord,lancuch
693 integer wartosc,default
696 iread=index(rekord,lancuch)
701 iread=iread+ilen(lancuch)+1
702 read (rekord(iread:),*) wartosc
705 C----------------------------------------------------------------------
706 subroutine multreadi(rekord,lancuch,tablica,dim,default)
709 integer tablica(dim),default
710 character*(*) rekord,lancuch
717 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
718 if (iread.eq.0) return
719 iread=iread+ilen(lancuch)+1
720 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
724 c----------------------------------------------------------------------------
725 subroutine card_concat(card)
727 include 'COMMON.IOUNITS'
729 character*80 karta,ucase
731 read (inp,'(a)') karta
734 do while (karta(80:80).eq.'&')
735 card=card(:ilen(card)+1)//karta(:79)
736 read (inp,'(a)') karta
739 card=card(:ilen(card)+1)//karta
742 c----------------------------------------------------------------------------
751 include 'COMMON.IOUNITS'
752 include 'COMMON.CONTROL'
753 integer lenpre,lenpot,ilen
755 character*16 cformat,cprint
757 integer lenint,lenout
758 call getenv('INPUT',prefix)
759 call getenv('OUTPUT',prefout)
760 call getenv('INTIN',prefintin)
761 call getenv('COORD',cformat)
762 call getenv('PRINTCOOR',cprint)
763 call getenv('SCRATCHDIR',scratchdir)
766 if (index(ucase(cformat),'CX').gt.0) then
773 lenint=ilen(prefintin)
774 C Get the names and open the input files
775 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
777 write (liczba,'(bz,i3.3)') me
778 outname=prefout(:lenout)//'_clust.out_'//liczba
780 outname=prefout(:lenout)//'_clust.out'
783 intinname=prefintin(:lenint)//'.bx'
784 else if (from_cx) then
785 intinname=prefintin(:lenint)//'.cx'
787 intinname=prefintin(:lenint)//'.int'
789 rmsname=prefintin(:lenint)//'.rms'
790 open (jplot,file=prefout(:ilen(prefout))//'.tex',
792 open (jrms,file=rmsname,status='unknown')
793 open(iout,file=outname,status='unknown')
794 C Get parameter filenames and open the parameter files.
795 call getenv('BONDPAR',bondname)
796 open (ibond,file=bondname,status='old')
797 call getenv('THETPAR',thetname)
798 open (ithep,file=thetname,status='old')
799 call getenv('ROTPAR',rotname)
800 open (irotam,file=rotname,status='old')
801 call getenv('TORPAR',torname)
802 open (itorp,file=torname,status='old')
803 call getenv('TORDPAR',tordname)
804 open (itordp,file=tordname,status='old')
805 call getenv('FOURIER',fouriername)
806 open (ifourier,file=fouriername,status='old')
807 call getenv('ELEPAR',elename)
808 open (ielep,file=elename,status='old')
809 call getenv('SIDEPAR',sidename)
810 open (isidep,file=sidename,status='old')
811 call getenv('SIDEP',sidepname)
812 open (isidep1,file=sidepname,status="old")
813 call getenv('SCCORPAR',sccorname)
814 open (isccor,file=sccorname,status="old")
815 call getenv('LIPTRANPAR',liptranname)
816 open (iliptranpar,file=liptranname,status='old')
819 C 8/9/01 In the newest version SCp interaction constants are read from a file
820 C Use -DOLDSCP to use hard-coded constants instead.
822 call getenv('SCPPAR',scpname)
823 open (iscpp,file=scpname,status='old')
827 subroutine read_dist_constr
828 implicit real*8 (a-h,o-z)
833 include 'COMMON.CONTROL'
834 include 'COMMON.CHAIN'
835 include 'COMMON.IOUNITS'
836 include 'COMMON.SBRIDGE'
837 integer ifrag_(2,100),ipair_(2,100)
838 double precision wfrag_(100),wpair_(100)
839 character*500 controlcard
840 logical lprn /.true./
841 write (iout,*) "Calling read_dist_constr"
842 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
844 write(iout,*) "TU sie wywalam?"
845 call card_concat(controlcard)
846 write (iout,*) controlcard
848 call readi(controlcard,"NFRAG",nfrag_,0)
849 call readi(controlcard,"NPAIR",npair_,0)
850 call readi(controlcard,"NDIST",ndist_,0)
851 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
852 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
853 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
854 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
855 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
856 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
857 write (iout,*) "IFRAG"
859 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
861 write (iout,*) "IPAIR"
863 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
867 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
868 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
869 & ifrag_(2,i)=nstart_sup+nsup-1
870 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
872 if (wfrag_(i).gt.0.0d0) then
873 do j=ifrag_(1,i),ifrag_(2,i)-1
875 write (iout,*) "j",j," k",k
877 if (constr_dist.eq.1) then
882 forcon(nhpb)=wfrag_(i)
883 else if (constr_dist.eq.2) then
884 if (ddjk.le.dist_cut) then
889 forcon(nhpb)=wfrag_(i)
896 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
899 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
900 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
906 if (wpair_(i).gt.0.0d0) then
914 do j=ifrag_(1,ii),ifrag_(2,ii)
915 do k=ifrag_(1,jj),ifrag_(2,jj)
919 forcon(nhpb)=wpair_(i)
921 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
922 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
928 if (constr_dist.eq.11) then
929 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
930 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
931 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
932 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
933 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
935 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
937 if (forcon(nhpb+1).gt.0.0d0) then
939 if (ibecarb(i).gt.0) then
943 if (dhpb(nhpb).eq.0.0d0)
944 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
945 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
946 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
947 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
960 c====-------------------------------------------------------------------
961 subroutine read_constr_homology(lprn)
967 include 'COMMON.SETUP'
968 include 'COMMON.CONTROL'
969 include 'COMMON.CHAIN'
970 include 'COMMON.IOUNITS'
972 include 'COMMON.INTERACT'
973 include 'COMMON.NAMES'
974 include 'COMMON.HOMRESTR'
979 c include 'include_unres/COMMON.VAR'
982 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
984 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
985 c & sigma_odl_temp(maxres,maxres,max_template)
987 character*24 model_ki_dist, model_ki_angle
988 character*500 controlcard
989 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
990 integer idomain(max_template,maxres)
994 logical unres_pdb,liiflag
996 c FP - Nov. 2014 Temporary specifications for new vars
998 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
1000 double precision, dimension (max_template,maxres) :: rescore
1001 double precision, dimension (max_template,maxres) :: rescore2
1002 double precision, dimension (max_template,maxres) :: rescore3
1003 character*24 tpl_k_rescore
1004 c -----------------------------------------------------------------
1005 c Reading multiple PDB ref structures and calculation of retraints
1006 c not using pre-computed ones stored in files model_ki_{dist,angle}
1008 c -----------------------------------------------------------------
1011 c Alternative: reading from input
1013 write (iout,*) "BEGIN READ HOMOLOGY INFO"
1020 call card_concat(controlcard)
1021 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1022 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1023 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1024 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1025 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1026 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1027 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1028 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1029 if (homol_nset.gt.1)then
1030 call readi(controlcard,"ISET",iset,1)
1031 call card_concat(controlcard)
1032 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1035 waga_homology(1)=1.0
1039 write(iout,*) "read_constr_homology iset",iset
1040 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1047 cd write (iout,*) "nnt",nnt," nct",nct
1057 c Reading HM global scores (prob not required)
1060 do k=1,constr_homology
1064 c open (4,file="HMscore")
1065 c do k=1,constr_homology
1066 c read (4,*,end=521) hmscore_tmp
1067 c hmscore(k)=hmscore_tmp ! Another transformation can be used
1068 c write(*,*) "Model", k, ":", hmscore(k)
1079 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1081 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1082 do k=1,constr_homology
1084 read(inp,'(a)') pdbfile
1085 c write (iout,*) "k ",k," pdbfile ",pdbfile
1086 c Next stament causes error upon compilation (?)
1087 c if(me.eq.king.or. .not. out1file)
1088 c write (iout,'(2a)') 'PDB data will be read from file ',
1089 c & pdbfile(:ilen(pdbfile))
1090 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1091 & pdbfile(:ilen(pdbfile))
1092 open(ipdbin,file=pdbfile,status='old',err=33)
1094 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1095 & pdbfile(:ilen(pdbfile))
1098 c print *,'Begin reading pdb data'
1100 c Files containing res sim or local scores (former containing sigmas)
1103 write(kic2,'(bz,i2.2)') k
1105 tpl_k_rescore="template"//kic2//".sco"
1108 call readpdb_template(k)
1111 crefjlee(j,i)=c(j,i)
1116 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1117 & (crefjlee(j,i+nres),j=1,3)
1119 write (iout,*) "READ HOMOLOGY INFO"
1120 write (iout,*) "read_constr_homology x: after reading pdb file"
1121 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1122 write (iout,*) "waga_dist",waga_dist
1123 write (iout,*) "waga_angle",waga_angle
1124 write (iout,*) "waga_theta",waga_theta
1125 write (iout,*) "waga_d",waga_d
1126 write (iout,*) "dist_cut",dist_cut
1135 c Distance restraints
1138 C Copy the coordinates from reference coordinates (?)
1142 c write (iout,*) "c(",j,i,") =",c(j,i)
1146 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1148 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1149 open (ientin,file=tpl_k_rescore,status='old')
1150 if (nnt.gt.1) rescore(k,1)=0.0d0
1151 do irec=nnt,nct ! loop for reading res sim
1152 if (read2sigma) then
1153 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1154 & rescore3_tmp,idomain_tmp
1156 idomain(k,i_tmp)=idomain_tmp
1157 rescore(k,i_tmp)=rescore_tmp
1158 rescore2(k,i_tmp)=rescore2_tmp
1159 rescore3(k,i_tmp)=rescore3_tmp
1160 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1161 & i_tmp,rescore2_tmp,rescore_tmp,
1162 & rescore3_tmp,idomain_tmp
1165 read (ientin,*,end=1401) rescore_tmp
1167 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1168 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1169 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1174 if (waga_dist.ne.0.0d0) then
1182 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1183 c write (iout,*) k,i,j,distal,dist2_cut
1185 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1186 & .and. distal.le.dist2_cut ) then
1192 c write (iout,*) "k",k
1193 c write (iout,*) "i",i," j",j," constr_homology",
1198 if (read2sigma) then
1201 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1203 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1204 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1205 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1207 if (odl(k,ii).le.dist_cut) then
1208 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1211 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1212 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1214 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1215 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1219 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1222 l_homo(k,ii)=.false.
1229 c Theta, dihedral and SC retraints
1231 if (waga_angle.gt.0.0d0) then
1232 c open (ientin,file=tpl_k_sigma_dih,status='old')
1233 c do irec=1,maxres-3 ! loop for reading sigma_dih
1234 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1235 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1236 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1237 c & sigma_dih(k,i+nnt-1)
1242 if (idomain(k,i).eq.0) then
1246 dih(k,i)=phiref(i) ! right?
1247 c read (ientin,*) sigma_dih(k,i) ! original variant
1248 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1249 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1250 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1251 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1252 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1254 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1255 & rescore(k,i-2)+rescore(k,i-3))/4.0
1256 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1257 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1258 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1259 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1260 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1261 c Instead of res sim other local measure of b/b str reliability possible
1262 if (sigma_dih(k,i).ne.0)
1263 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1264 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1269 if (waga_theta.gt.0.0d0) then
1270 c open (ientin,file=tpl_k_sigma_theta,status='old')
1271 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1272 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1273 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1274 c & sigma_theta(k,i+nnt-1)
1279 do i = nnt+2,nct ! right? without parallel.
1280 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1281 c do i=ithet_start,ithet_end ! with FG parallel.
1282 if (idomain(k,i).eq.0) then
1283 sigma_theta(k,i)=0.0
1286 thetatpl(k,i)=thetaref(i)
1287 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1288 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1289 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1290 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1291 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1292 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1293 & rescore(k,i-2))/3.0
1294 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1295 if (sigma_theta(k,i).ne.0)
1296 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1298 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1299 c rescore(k,i-2) ! right expression ?
1300 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1304 if (waga_d.gt.0.0d0) then
1305 c open (ientin,file=tpl_k_sigma_d,status='old')
1306 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1307 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1308 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1309 c & sigma_d(k,i+nnt-1)
1313 do i = nnt,nct ! right? without parallel.
1314 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1315 c do i=loc_start,loc_end ! with FG parallel.
1316 if (itype(i).eq.10) cycle
1317 if (idomain(k,i).eq.0 ) then
1324 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1325 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1326 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1327 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1328 sigma_d(k,i)=rescore3(k,i) ! right expression ?
1329 if (sigma_d(k,i).ne.0)
1330 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1332 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1333 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1334 c read (ientin,*) sigma_d(k,i) ! 1st variant
1339 c remove distance restraints not used in any model from the list
1340 c shift data in all arrays
1342 if (waga_dist.ne.0.0d0) then
1348 if (ii_in_use(ii).eq.0.and.liiflag) then
1352 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1353 & .not.liiflag.and.ii.eq.lim_odl) then
1354 if (ii.eq.lim_odl) then
1355 iishift=ii-iistart+1
1360 do ki=iistart,lim_odl-iishift
1361 ires_homo(ki)=ires_homo(ki+iishift)
1362 jres_homo(ki)=jres_homo(ki+iishift)
1363 ii_in_use(ki)=ii_in_use(ki+iishift)
1364 do k=1,constr_homology
1365 odl(k,ki)=odl(k,ki+iishift)
1366 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1367 l_homo(k,ki)=l_homo(k,ki+iishift)
1371 lim_odl=lim_odl-iishift
1376 if (constr_homology.gt.0) call homology_partition
1377 if (constr_homology.gt.0) call init_int_table
1378 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1379 cd & "lim_xx=",lim_xx
1380 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1381 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1385 if (.not.lprn) return
1386 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1387 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1388 write (iout,*) "Distance restraints from templates"
1390 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1391 & ii,ires_homo(ii),jres_homo(ii),
1392 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1393 & ki=1,constr_homology)
1395 write (iout,*) "Dihedral angle restraints from templates"
1397 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1398 & (rad2deg*dih(ki,i),
1399 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1401 write (iout,*) "Virtual-bond angle restraints from templates"
1403 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1404 & (rad2deg*thetatpl(ki,i),
1405 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1407 write (iout,*) "SC restraints from templates"
1409 write(iout,'(i5,100(4f8.2,4x))') i,
1410 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1411 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1414 c -----------------------------------------------------------------