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 call reada(weightcard,'WSC',wsc,1.0d0)
172 call reada(weightcard,'WLONG',wsc,wsc)
173 call reada(weightcard,'WSCP',wscp,1.0d0)
174 call reada(weightcard,'WELEC',welec,1.0D0)
175 call reada(weightcard,'WVDWPP',wvdwpp,welec)
176 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
177 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
178 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
179 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
180 call reada(weightcard,'WTURN3',wturn3,1.0D0)
181 call reada(weightcard,'WTURN4',wturn4,1.0D0)
182 call reada(weightcard,'WTURN6',wturn6,1.0D0)
183 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
184 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
185 call reada(weightcard,'WBOND',wbond,1.0D0)
186 call reada(weightcard,'WTOR',wtor,1.0D0)
187 call reada(weightcard,'WTORD',wtor_d,1.0D0)
188 call reada(weightcard,'WANG',wang,1.0D0)
189 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
190 call reada(weightcard,'SCAL14',scal14,0.4D0)
191 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
192 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
193 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
194 if (index(weightcard,'SOFT').gt.0) ipot=6
195 call reada(weightcard,"D0CM",d0cm,3.78d0)
196 call reada(weightcard,"AKCM",akcm,15.1d0)
197 call reada(weightcard,"AKTH",akth,11.0d0)
198 call reada(weightcard,"AKCT",akct,12.0d0)
199 call reada(weightcard,"V1SS",v1ss,-1.08d0)
200 call reada(weightcard,"V2SS",v2ss,7.61d0)
201 call reada(weightcard,"V3SS",v3ss,13.7d0)
202 call reada(weightcard,"EBR",ebr,-5.50D0)
203 call reada(weightcard,'WSHIELD',wshield,1.0d0)
204 write(iout,*) 'WSHIELD',wshield
205 call reada(weightcard,'WLT',wliptran,0.0D0)
206 call reada(weightcard,"ATRISS",atriss,0.301D0)
207 call reada(weightcard,"BTRISS",btriss,0.021D0)
208 call reada(weightcard,"CTRISS",ctriss,1.001D0)
209 call reada(weightcard,"DTRISS",dtriss,1.001D0)
210 write (iout,*) "ATRISS=", atriss
211 write (iout,*) "BTRISS=", btriss
212 write (iout,*) "CTRISS=", ctriss
213 write (iout,*) "DTRISS=", dtriss
214 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
216 dyn_ss_mask(i)=.false.
220 dyn_ssbond_ij(i,j)=1.0d300
223 call reada(weightcard,"HT",Ht,0.0D0)
225 ss_depth=ebr/wsc-0.25*eps(1,1)
226 Ht=Ht/wsc-0.25*eps(1,1)
227 akcm=akcm*wstrain/wsc
228 akth=akth*wstrain/wsc
229 akct=akct*wstrain/wsc
230 v1ss=v1ss*wstrain/wsc
231 v2ss=v2ss*wstrain/wsc
232 v3ss=v3ss*wstrain/wsc
234 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
236 write (iout,'(/a)') "Disulfide bridge parameters:"
237 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
238 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
239 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
240 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
243 C 12/1/95 Added weight for the multi-body term WCORR
244 call reada(weightcard,'WCORRH',wcorr,1.0D0)
245 if (wcorr4.gt.0.0d0) wcorr=wcorr4
265 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
266 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
267 & wturn4,wturn6,wsccor
268 10 format (/'Energy-term weights (unscaled):'//
269 & 'WSCC= ',f10.6,' (SC-SC)'/
270 & 'WSCP= ',f10.6,' (SC-p)'/
271 & 'WELEC= ',f10.6,' (p-p electr)'/
272 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
273 & 'WBOND= ',f10.6,' (stretching)'/
274 & 'WANG= ',f10.6,' (bending)'/
275 & 'WSCLOC= ',f10.6,' (SC local)'/
276 & 'WTOR= ',f10.6,' (torsional)'/
277 & 'WTORD= ',f10.6,' (double torsional)'/
278 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
279 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
280 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
281 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
282 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
283 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
284 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
285 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
286 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
288 if (wcorr4.gt.0.0d0) then
289 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
290 & 'between contact pairs of peptide groups'
291 write (iout,'(2(a,f5.3/))')
292 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
293 & 'Range of quenching the correlation terms:',2*delt_corr
294 else if (wcorr.gt.0.0d0) then
295 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
296 & 'between contact pairs of peptide groups'
298 write (iout,'(a,f8.3)')
299 & 'Scaling factor of 1,4 SC-p interactions:',scal14
300 write (iout,'(a,f8.3)')
301 & 'General scaling factor of SC-p interactions:',scalscp
302 r0_corr=cutoff_corr-delt_corr
304 aad(i,1)=scalscp*aad(i,1)
305 aad(i,2)=scalscp*aad(i,2)
306 bad(i,1)=scalscp*bad(i,1)
307 bad(i,2)=scalscp*bad(i,2)
311 print *,'indpdb=',indpdb,' pdbref=',pdbref
313 C Read sequence if not taken from the pdb file.
314 if (iscode.gt.0) then
315 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
317 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
319 C Convert sequence to numeric code
321 itype(i)=rescode(i,sequence(i),iscode)
324 print '(20i4)',(itype(i),i=1,nres)
328 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
330 if (itype(i).eq.ntyp1) then
334 else if (iabs(itype(i+1)).ne.20) then
336 else if (iabs(itype(i)).ne.20) then
343 write (iout,*) "ITEL"
345 write (iout,*) i,itype(i),itel(i)
348 print *,'Call Read_Bridge.'
350 C this fragment reads diheadral constrains
351 if (with_dihed_constr) then
353 read (inp,*) ndih_constr
354 if (ndih_constr.gt.0) then
356 C write (iout,*) 'FTORS',ftors
357 C ftors is the force constant for torsional quartic constrains
358 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
361 & 'There are',ndih_constr,' constraints on phi angles.'
363 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
367 phi0(i)=deg2rad*phi0(i)
368 drange(i)=deg2rad*drange(i)
370 endif ! endif ndif_constr.gt.0
371 endif ! with_dihed_constr
372 if (with_theta_constr) then
373 C with_theta_constr is keyword allowing for occurance of theta constrains
374 read (inp,*) ntheta_constr
375 C ntheta_constr is the number of theta constrains
376 if (ntheta_constr.gt.0) then
378 read (inp,*) (itheta_constr(i),theta_constr0(i),
379 & theta_drange(i),for_thet_constr(i),
381 C the above code reads from 1 to ntheta_constr
382 C itheta_constr(i) residue i for which is theta_constr
383 C theta_constr0 the global minimum value
384 C theta_drange is range for which there is no energy penalty
385 C for_thet_constr is the force constant for quartic energy penalty
387 C if(me.eq.king.or..not.out1file)then
389 & 'There are',ntheta_constr,' constraints on phi angles.'
391 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
397 theta_constr0(i)=deg2rad*theta_constr0(i)
398 theta_drange(i)=deg2rad*theta_drange(i)
400 C if(me.eq.king.or..not.out1file)
401 C & write (iout,*) 'FTORS',ftors
402 C do i=1,ntheta_constr
403 C ii = itheta_constr(i)
404 C thetabound(1,ii) = phi0(i)-drange(i)
405 C thetabound(2,ii) = phi0(i)+drange(i)
407 endif ! ntheta_constr.gt.0
408 endif! with_theta_constr
412 print *,'NNT=',NNT,' NCT=',NCT
413 if (itype(1).eq.ntyp1) nnt=2
414 if (itype(nres).eq.ntyp1) nct=nct-1
415 if (nstart.lt.nnt) nstart=nnt
416 if (nend.gt.nct .or. nend.eq.0) nend=nct
417 write (iout,*) "nstart",nstart," nend",nend
419 if (constr_homology.gt.0) then
420 call read_constr_homology(print_homology_restraints)
424 c read(inp,'(a)') pdbfile
425 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
426 c open(ipdbin,file=pdbfile,status='old',err=33)
428 c 33 write (iout,'(a)') 'Error opening PDB file.'
431 c print *,'Begin reading pdb data'
433 c print *,'Finished reading pdb data'
434 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
436 c itype_pdb(i)=itype(i)
439 c write (iout,'(a,i3)') 'nsup=',nsup
441 c if (nsup.le.(nct-nnt+1)) then
442 c do i=0,nct-nnt+1-nsup
443 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
449 c & 'Error - sequences to be superposed do not match.'
452 c do i=0,nsup-(nct-nnt+1)
453 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
455 c nstart_sup=nstart_sup+i
461 c & 'Error - sequences to be superposed do not match.'
464 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
465 c & ' nstart_seq=',nstart_seq
469 write (iout,*) "molread: REFSTR",refstr
471 if (.not.pdbref) then
472 call read_angles(inp,*38)
474 38 write (iout,'(a)') 'Error reading reference structure.'
476 call mp_stopall(Error_Msg)
478 stop 'Error reading reference structure'
491 call contact(.true.,ncont_ref,icont_ref)
494 C write (iout,'(/a,i3,a)')
495 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
496 write (iout,'(20i4)') (iss(i),i=1,ns)
498 write(iout,*)"Running with dynamic disulfide-bond formation"
500 write (iout,'(/a/)') 'Pre-formed links are:'
506 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
507 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
513 if (ns.gt.0.and.dyn_ss) then
517 forcon(i-nss)=forcon(i)
524 dyn_ss_mask(iss(i))=.true.
527 c Read distance restraints
528 if (constr_dist.gt.0) then
529 call read_dist_constr
534 c-----------------------------------------------------------------------------
535 logical function seq_comp(itypea,itypeb,length)
537 integer length,itypea(length),itypeb(length)
540 if (itypea(i).ne.itypeb(i)) then
548 c-----------------------------------------------------------------------------
549 subroutine read_bridge
550 C Read information about disulfide bridges.
553 include 'COMMON.IOUNITS'
556 include 'COMMON.INTERACT'
557 include 'COMMON.LOCAL'
558 include 'COMMON.NAMES'
559 include 'COMMON.CHAIN'
560 include 'COMMON.FFIELD'
561 include 'COMMON.SBRIDGE'
562 include 'COMMON.HEADER'
563 include 'COMMON.CONTROL'
564 include 'COMMON.TIME1'
566 include 'COMMON.INFO'
569 C Read bridging residues.
570 read (inp,*) ns,(iss(i),i=1,ns)
572 C Check whether the specified bridging residues are cystines.
574 if (itype(iss(i)).ne.1) then
575 write (iout,'(2a,i3,a)')
576 & 'Do you REALLY think that the residue ',
577 & restyp(itype(iss(i))),i,
578 & ' can form a disulfide bridge?!!!'
579 write (*,'(2a,i3,a)')
580 & 'Do you REALLY think that the residue ',
581 & restyp(itype(iss(i))),i,
582 & ' can form a disulfide bridge?!!!'
584 call mp_stopall(error_msg)
590 C Read preformed bridges.
592 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
595 C Check if the residues involved in bridges are in the specified list of
599 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
600 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
601 write (iout,'(a,i3,a)') 'Disulfide pair',i,
602 & ' contains residues present in other pairs.'
603 write (*,'(a,i3,a)') 'Disulfide pair',i,
604 & ' contains residues present in other pairs.'
606 call mp_stopall(error_msg)
613 if (ihpb(i).eq.iss(j)) goto 10
615 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
618 if (jhpb(i).eq.iss(j)) goto 20
620 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
633 c----------------------------------------------------------------------------
634 subroutine read_angles(kanal,*)
639 include 'COMMON.CHAIN'
640 include 'COMMON.IOUNITS'
642 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
643 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
644 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
645 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
647 theta(i)=deg2rad*theta(i)
648 phi(i)=deg2rad*phi(i)
649 alph(i)=deg2rad*alph(i)
650 omeg(i)=deg2rad*omeg(i)
655 c----------------------------------------------------------------------------
656 subroutine reada(rekord,lancuch,wartosc,default)
658 character*(*) rekord,lancuch
659 double precision wartosc,default
662 iread=index(rekord,lancuch)
667 iread=iread+ilen(lancuch)+1
668 read (rekord(iread:),*) wartosc
671 c----------------------------------------------------------------------------
672 subroutine multreada(rekord,lancuch,tablica,dim,default)
675 double precision tablica(dim),default
676 character*(*) rekord,lancuch
682 iread=index(rekord,lancuch)
683 if (iread.eq.0) return
684 iread=iread+ilen(lancuch)+1
685 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
688 c----------------------------------------------------------------------------
689 subroutine readi(rekord,lancuch,wartosc,default)
691 character*(*) rekord,lancuch
692 integer wartosc,default
695 iread=index(rekord,lancuch)
700 iread=iread+ilen(lancuch)+1
701 read (rekord(iread:),*) wartosc
704 C----------------------------------------------------------------------
705 subroutine multreadi(rekord,lancuch,tablica,dim,default)
708 integer tablica(dim),default
709 character*(*) rekord,lancuch
716 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
717 if (iread.eq.0) return
718 iread=iread+ilen(lancuch)+1
719 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
723 c----------------------------------------------------------------------------
724 subroutine card_concat(card)
726 include 'COMMON.IOUNITS'
728 character*80 karta,ucase
730 read (inp,'(a)') karta
733 do while (karta(80:80).eq.'&')
734 card=card(:ilen(card)+1)//karta(:79)
735 read (inp,'(a)') karta
738 card=card(:ilen(card)+1)//karta
741 c----------------------------------------------------------------------------
750 include 'COMMON.IOUNITS'
751 include 'COMMON.CONTROL'
752 integer lenpre,lenpot,ilen
754 character*16 cformat,cprint
756 integer lenint,lenout
757 call getenv('INPUT',prefix)
758 call getenv('OUTPUT',prefout)
759 call getenv('INTIN',prefintin)
760 call getenv('COORD',cformat)
761 call getenv('PRINTCOOR',cprint)
762 call getenv('SCRATCHDIR',scratchdir)
765 if (index(ucase(cformat),'CX').gt.0) then
772 lenint=ilen(prefintin)
773 C Get the names and open the input files
774 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
776 write (liczba,'(bz,i3.3)') me
777 outname=prefout(:lenout)//'_clust.out_'//liczba
779 outname=prefout(:lenout)//'_clust.out'
782 intinname=prefintin(:lenint)//'.bx'
783 else if (from_cx) then
784 intinname=prefintin(:lenint)//'.cx'
786 intinname=prefintin(:lenint)//'.int'
788 rmsname=prefintin(:lenint)//'.rms'
789 open (jplot,file=prefout(:ilen(prefout))//'.tex',
791 open (jrms,file=rmsname,status='unknown')
792 open(iout,file=outname,status='unknown')
793 C Get parameter filenames and open the parameter files.
794 call getenv('BONDPAR',bondname)
795 open (ibond,file=bondname,status='old')
796 call getenv('THETPAR',thetname)
797 open (ithep,file=thetname,status='old')
798 call getenv('ROTPAR',rotname)
799 open (irotam,file=rotname,status='old')
800 call getenv('TORPAR',torname)
801 open (itorp,file=torname,status='old')
802 call getenv('TORDPAR',tordname)
803 open (itordp,file=tordname,status='old')
804 call getenv('FOURIER',fouriername)
805 open (ifourier,file=fouriername,status='old')
806 call getenv('ELEPAR',elename)
807 open (ielep,file=elename,status='old')
808 call getenv('SIDEPAR',sidename)
809 open (isidep,file=sidename,status='old')
810 call getenv('SIDEP',sidepname)
811 open (isidep1,file=sidepname,status="old")
812 call getenv('SCCORPAR',sccorname)
813 open (isccor,file=sccorname,status="old")
814 call getenv('LIPTRANPAR',liptranname)
815 open (iliptranpar,file=liptranname,status='old')
818 C 8/9/01 In the newest version SCp interaction constants are read from a file
819 C Use -DOLDSCP to use hard-coded constants instead.
821 call getenv('SCPPAR',scpname)
822 open (iscpp,file=scpname,status='old')
826 subroutine read_dist_constr
827 implicit real*8 (a-h,o-z)
832 include 'COMMON.CONTROL'
833 include 'COMMON.CHAIN'
834 include 'COMMON.IOUNITS'
835 include 'COMMON.SBRIDGE'
836 integer ifrag_(2,100),ipair_(2,100)
837 double precision wfrag_(100),wpair_(100)
838 character*500 controlcard
839 logical lprn /.true./
840 write (iout,*) "Calling read_dist_constr"
841 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
843 write(iout,*) "TU sie wywalam?"
844 call card_concat(controlcard)
845 write (iout,*) controlcard
847 call readi(controlcard,"NFRAG",nfrag_,0)
848 call readi(controlcard,"NPAIR",npair_,0)
849 call readi(controlcard,"NDIST",ndist_,0)
850 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
851 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
852 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
853 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
854 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
855 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
856 write (iout,*) "IFRAG"
858 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
860 write (iout,*) "IPAIR"
862 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
866 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
867 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
868 & ifrag_(2,i)=nstart_sup+nsup-1
869 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
871 if (wfrag_(i).gt.0.0d0) then
872 do j=ifrag_(1,i),ifrag_(2,i)-1
874 write (iout,*) "j",j," k",k
876 if (constr_dist.eq.1) then
881 forcon(nhpb)=wfrag_(i)
882 else if (constr_dist.eq.2) then
883 if (ddjk.le.dist_cut) then
888 forcon(nhpb)=wfrag_(i)
895 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
898 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
899 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
905 if (wpair_(i).gt.0.0d0) then
913 do j=ifrag_(1,ii),ifrag_(2,ii)
914 do k=ifrag_(1,jj),ifrag_(2,jj)
918 forcon(nhpb)=wpair_(i)
920 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
921 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
927 if (constr_dist.eq.11) then
928 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
929 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
930 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
931 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
932 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
934 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
936 if (forcon(nhpb+1).gt.0.0d0) then
938 if (ibecarb(i).gt.0) then
942 if (dhpb(nhpb).eq.0.0d0)
943 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
944 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
945 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
946 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
959 c====-------------------------------------------------------------------
960 subroutine read_constr_homology(lprn)
966 include 'COMMON.SETUP'
967 include 'COMMON.CONTROL'
968 include 'COMMON.CHAIN'
969 include 'COMMON.IOUNITS'
971 include 'COMMON.INTERACT'
972 include 'COMMON.NAMES'
973 include 'COMMON.HOMRESTR'
978 c include 'include_unres/COMMON.VAR'
981 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
983 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
984 c & sigma_odl_temp(maxres,maxres,max_template)
986 character*24 model_ki_dist, model_ki_angle
987 character*500 controlcard
988 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
989 integer idomain(max_template,maxres)
993 logical unres_pdb,liiflag
995 c FP - Nov. 2014 Temporary specifications for new vars
997 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
999 double precision, dimension (max_template,maxres) :: rescore
1000 double precision, dimension (max_template,maxres) :: rescore2
1001 double precision, dimension (max_template,maxres) :: rescore3
1002 character*24 tpl_k_rescore
1003 c -----------------------------------------------------------------
1004 c Reading multiple PDB ref structures and calculation of retraints
1005 c not using pre-computed ones stored in files model_ki_{dist,angle}
1007 c -----------------------------------------------------------------
1010 c Alternative: reading from input
1012 write (iout,*) "BEGIN READ HOMOLOGY INFO"
1019 call card_concat(controlcard)
1020 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1021 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1022 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1023 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1024 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1025 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1026 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1027 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1028 if (homol_nset.gt.1)then
1029 call readi(controlcard,"ISET",iset,1)
1030 call card_concat(controlcard)
1031 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1034 waga_homology(1)=1.0
1038 write(iout,*) "read_constr_homology iset",iset
1039 write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1046 cd write (iout,*) "nnt",nnt," nct",nct
1056 c Reading HM global scores (prob not required)
1059 do k=1,constr_homology
1063 c open (4,file="HMscore")
1064 c do k=1,constr_homology
1065 c read (4,*,end=521) hmscore_tmp
1066 c hmscore(k)=hmscore_tmp ! Another transformation can be used
1067 c write(*,*) "Model", k, ":", hmscore(k)
1078 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1080 write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1081 do k=1,constr_homology
1083 read(inp,'(a)') pdbfile
1084 c write (iout,*) "k ",k," pdbfile ",pdbfile
1085 c Next stament causes error upon compilation (?)
1086 c if(me.eq.king.or. .not. out1file)
1087 c write (iout,'(2a)') 'PDB data will be read from file ',
1088 c & pdbfile(:ilen(pdbfile))
1089 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1090 & pdbfile(:ilen(pdbfile))
1091 open(ipdbin,file=pdbfile,status='old',err=33)
1093 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1094 & pdbfile(:ilen(pdbfile))
1097 c print *,'Begin reading pdb data'
1099 c Files containing res sim or local scores (former containing sigmas)
1102 write(kic2,'(bz,i2.2)') k
1104 tpl_k_rescore="template"//kic2//".sco"
1107 call readpdb_template(k)
1110 crefjlee(j,i)=c(j,i)
1115 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1116 & (crefjlee(j,i+nres),j=1,3)
1118 write (iout,*) "READ HOMOLOGY INFO"
1119 write (iout,*) "read_constr_homology x: after reading pdb file"
1120 write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1121 write (iout,*) "waga_dist",waga_dist
1122 write (iout,*) "waga_angle",waga_angle
1123 write (iout,*) "waga_theta",waga_theta
1124 write (iout,*) "waga_d",waga_d
1125 write (iout,*) "dist_cut",dist_cut
1134 c Distance restraints
1137 C Copy the coordinates from reference coordinates (?)
1141 c write (iout,*) "c(",j,i,") =",c(j,i)
1145 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1147 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1148 open (ientin,file=tpl_k_rescore,status='old')
1149 if (nnt.gt.1) rescore(k,1)=0.0d0
1150 do irec=nnt,nct ! loop for reading res sim
1151 if (read2sigma) then
1152 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1153 & rescore3_tmp,idomain_tmp
1155 idomain(k,i_tmp)=idomain_tmp
1156 rescore(k,i_tmp)=rescore_tmp
1157 rescore2(k,i_tmp)=rescore2_tmp
1158 rescore3(k,i_tmp)=rescore3_tmp
1159 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1160 & i_tmp,rescore2_tmp,rescore_tmp,
1161 & rescore3_tmp,idomain_tmp
1164 read (ientin,*,end=1401) rescore_tmp
1166 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1167 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1168 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1173 if (waga_dist.ne.0.0d0) then
1181 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1182 c write (iout,*) k,i,j,distal,dist2_cut
1184 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1185 & .and. distal.le.dist2_cut ) then
1191 c write (iout,*) "k",k
1192 c write (iout,*) "i",i," j",j," constr_homology",
1197 if (read2sigma) then
1200 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1202 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1203 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1204 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1206 if (odl(k,ii).le.dist_cut) then
1207 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1210 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1211 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1213 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1214 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1218 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1221 l_homo(k,ii)=.false.
1228 c Theta, dihedral and SC retraints
1230 if (waga_angle.gt.0.0d0) then
1231 c open (ientin,file=tpl_k_sigma_dih,status='old')
1232 c do irec=1,maxres-3 ! loop for reading sigma_dih
1233 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1234 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1235 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1236 c & sigma_dih(k,i+nnt-1)
1241 if (idomain(k,i).eq.0) then
1245 dih(k,i)=phiref(i) ! right?
1246 c read (ientin,*) sigma_dih(k,i) ! original variant
1247 c write (iout,*) "dih(",k,i,") =",dih(k,i)
1248 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1249 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1250 c & "rescore(",k,i-2,") =",rescore(k,i-2),
1251 c & "rescore(",k,i-3,") =",rescore(k,i-3)
1253 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1254 & rescore(k,i-2)+rescore(k,i-3))/4.0
1255 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1256 c write (iout,*) "Raw sigmas for dihedral angle restraints"
1257 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1258 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1259 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1260 c Instead of res sim other local measure of b/b str reliability possible
1261 if (sigma_dih(k,i).ne.0)
1262 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1263 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1268 if (waga_theta.gt.0.0d0) then
1269 c open (ientin,file=tpl_k_sigma_theta,status='old')
1270 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1271 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1272 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1273 c & sigma_theta(k,i+nnt-1)
1278 do i = nnt+2,nct ! right? without parallel.
1279 c do i = i=1,nres ! alternative for bounds acc to readpdb?
1280 c do i=ithet_start,ithet_end ! with FG parallel.
1281 if (idomain(k,i).eq.0) then
1282 sigma_theta(k,i)=0.0
1285 thetatpl(k,i)=thetaref(i)
1286 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1287 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1288 c & "rescore(",k,i-1,") =",rescore(k,i-1),
1289 c & "rescore(",k,i-2,") =",rescore(k,i-2)
1290 c read (ientin,*) sigma_theta(k,i) ! 1st variant
1291 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1292 & rescore(k,i-2))/3.0
1293 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1294 if (sigma_theta(k,i).ne.0)
1295 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1297 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1298 c rescore(k,i-2) ! right expression ?
1299 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1303 if (waga_d.gt.0.0d0) then
1304 c open (ientin,file=tpl_k_sigma_d,status='old')
1305 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1306 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1307 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1308 c & sigma_d(k,i+nnt-1)
1312 do i = nnt,nct ! right? without parallel.
1313 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
1314 c do i=loc_start,loc_end ! with FG parallel.
1315 if (itype(i).eq.10) cycle
1316 if (idomain(k,i).eq.0 ) then
1323 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1324 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1325 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1326 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1327 sigma_d(k,i)=rescore3(k,i) ! right expression ?
1328 if (sigma_d(k,i).ne.0)
1329 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1331 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1332 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1333 c read (ientin,*) sigma_d(k,i) ! 1st variant
1338 c remove distance restraints not used in any model from the list
1339 c shift data in all arrays
1341 if (waga_dist.ne.0.0d0) then
1347 if (ii_in_use(ii).eq.0.and.liiflag) then
1351 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1352 & .not.liiflag.and.ii.eq.lim_odl) then
1353 if (ii.eq.lim_odl) then
1354 iishift=ii-iistart+1
1359 do ki=iistart,lim_odl-iishift
1360 ires_homo(ki)=ires_homo(ki+iishift)
1361 jres_homo(ki)=jres_homo(ki+iishift)
1362 ii_in_use(ki)=ii_in_use(ki+iishift)
1363 do k=1,constr_homology
1364 odl(k,ki)=odl(k,ki+iishift)
1365 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1366 l_homo(k,ki)=l_homo(k,ki+iishift)
1370 lim_odl=lim_odl-iishift
1375 if (constr_homology.gt.0) call homology_partition
1376 if (constr_homology.gt.0) call init_int_table
1377 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1378 cd & "lim_xx=",lim_xx
1379 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1380 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1384 if (.not.lprn) return
1385 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1386 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1387 write (iout,*) "Distance restraints from templates"
1389 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1390 & ii,ires_homo(ii),jres_homo(ii),
1391 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1392 & ki=1,constr_homology)
1394 write (iout,*) "Dihedral angle restraints from templates"
1396 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1397 & (rad2deg*dih(ki,i),
1398 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1400 write (iout,*) "Virtual-bond angle restraints from templates"
1402 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1403 & (rad2deg*thetatpl(ki,i),
1404 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1406 write (iout,*) "SC restraints from templates"
1408 write(iout,'(i5,100(4f8.2,4x))') i,
1409 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1410 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1413 c -----------------------------------------------------------------