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 character*320 controlcard,ucase
22 integer i,i1,i2,it1,it2
24 read (INP,'(a80)') titel
25 call card_concat(controlcard)
27 call readi(controlcard,'NRES',nres,0)
28 call readi(controlcard,'RESCALE',rescale_mode,2)
29 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
30 write (iout,*) "DISTCHAINMAX",distchainmax
31 call readi(controlcard,'PDBOUT',outpdb,0)
32 call readi(controlcard,'MOL2OUT',outmol2,0)
33 refstr=(index(controlcard,'REFSTR').gt.0)
34 write (iout,*) "REFSTR",refstr
35 pdbref=(index(controlcard,'PDBREF').gt.0)
36 iscode=index(controlcard,'ONE_LETTER')
37 tree=(index(controlcard,'MAKE_TREE').gt.0)
38 min_var=(index(controlcard,'MINVAR').gt.0)
39 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
40 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
41 call readi(controlcard,'NCUT',ncut,1)
42 call readi(controlcard,'SYM',symetr,1)
43 write (iout,*) 'sym', symetr
44 call readi(controlcard,'NSTART',nstart,0)
45 call readi(controlcard,'NEND',nend,0)
46 call reada(controlcard,'ECUT',ecut,10.0d0)
47 call reada(controlcard,'PROB',prob_limit,0.99d0)
48 write (iout,*) "Probability limit",prob_limit
49 lgrp=(index(controlcard,'LGRP').gt.0)
50 caonly=(index(controlcard,'CA_ONLY').gt.0)
51 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
52 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
53 call readi(controlcard,'IOPT',iopt,2)
54 lside = index(controlcard,"SIDE").gt.0
55 efree = index(controlcard,"EFREE").gt.0
56 call readi(controlcard,'NTEMP',nT,1)
57 write (iout,*) "nT",nT
58 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
59 write (iout,*) "nT",nT
60 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
62 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
64 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
65 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
66 lprint_int=index(controlcard,"PRINT_INT") .gt.0
70 c--------------------------------------------------------------------------
73 C Read molecular data.
77 include 'COMMON.IOUNITS'
80 include 'COMMON.INTERACT'
81 include 'COMMON.LOCAL'
82 include 'COMMON.NAMES'
83 include 'COMMON.CHAIN'
84 include 'COMMON.FFIELD'
85 include 'COMMON.SBRIDGE'
86 include 'COMMON.HEADER'
87 include 'COMMON.CONTROL'
88 include 'COMMON.CONTACTS'
89 include 'COMMON.TIME1'
93 character*4 sequence(maxres)
94 character*800 weightcard
96 double precision x(maxvar)
97 integer itype_pdb(maxres)
99 integer i,j,kkk,i1,i2,it1,it2
103 C Read weights of the subsequent energy terms.
104 call card_concat(weightcard)
105 call reada(weightcard,'WSC',wsc,1.0d0)
106 call reada(weightcard,'WLONG',wsc,wsc)
107 call reada(weightcard,'WSCP',wscp,1.0d0)
108 call reada(weightcard,'WELEC',welec,1.0D0)
109 call reada(weightcard,'WVDWPP',wvdwpp,welec)
110 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
111 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
112 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
113 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
114 call reada(weightcard,'WTURN3',wturn3,1.0D0)
115 call reada(weightcard,'WTURN4',wturn4,1.0D0)
116 call reada(weightcard,'WTURN6',wturn6,1.0D0)
117 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
118 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
119 call reada(weightcard,'WBOND',wbond,1.0D0)
120 call reada(weightcard,'WTOR',wtor,1.0D0)
121 call reada(weightcard,'WTORD',wtor_d,1.0D0)
122 call reada(weightcard,'WANG',wang,1.0D0)
123 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
124 call reada(weightcard,'SCAL14',scal14,0.4D0)
125 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
126 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
127 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
128 if (index(weightcard,'SOFT').gt.0) ipot=6
129 call reada(weightcard,"D0CM",d0cm,3.78d0)
130 call reada(weightcard,"AKCM",akcm,15.1d0)
131 call reada(weightcard,"AKTH",akth,11.0d0)
132 call reada(weightcard,"AKCT",akct,12.0d0)
133 call reada(weightcard,"V1SS",v1ss,-1.08d0)
134 call reada(weightcard,"V2SS",v2ss,7.61d0)
135 call reada(weightcard,"V3SS",v3ss,13.7d0)
136 call reada(weightcard,"EBR",ebr,-5.50D0)
137 call reada(weightcard,"ATRISS",atriss,0.301D0)
138 call reada(weightcard,"BTRISS",btriss,0.021D0)
139 call reada(weightcard,"CTRISS",ctriss,1.001D0)
140 call reada(weightcard,"DTRISS",dtriss,1.001D0)
141 write (iout,*) "ATRISS=", atriss
142 write (iout,*) "BTRISS=", btriss
143 write (iout,*) "CTRISS=", ctriss
144 write (iout,*) "DTRISS=", dtriss
145 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
147 dyn_ss_mask(i)=.false.
151 dyn_ssbond_ij(i,j)=1.0d300
154 call reada(weightcard,"HT",Ht,0.0D0)
156 ss_depth=ebr/wsc-0.25*eps(1,1)
157 Ht=Ht/wsc-0.25*eps(1,1)
158 akcm=akcm*wstrain/wsc
159 akth=akth*wstrain/wsc
160 akct=akct*wstrain/wsc
161 v1ss=v1ss*wstrain/wsc
162 v2ss=v2ss*wstrain/wsc
163 v3ss=v3ss*wstrain/wsc
165 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
167 write (iout,'(/a)') "Disulfide bridge parameters:"
168 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
169 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
170 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
171 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
174 C 12/1/95 Added weight for the multi-body term WCORR
175 call reada(weightcard,'WCORRH',wcorr,1.0D0)
176 if (wcorr4.gt.0.0d0) wcorr=wcorr4
196 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
197 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
198 & wturn4,wturn6,wsccor
199 10 format (/'Energy-term weights (unscaled):'//
200 & 'WSCC= ',f10.6,' (SC-SC)'/
201 & 'WSCP= ',f10.6,' (SC-p)'/
202 & 'WELEC= ',f10.6,' (p-p electr)'/
203 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
204 & 'WBOND= ',f10.6,' (stretching)'/
205 & 'WANG= ',f10.6,' (bending)'/
206 & 'WSCLOC= ',f10.6,' (SC local)'/
207 & 'WTOR= ',f10.6,' (torsional)'/
208 & 'WTORD= ',f10.6,' (double torsional)'/
209 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
210 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
211 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
212 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
213 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
214 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
215 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
216 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
217 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
219 if (wcorr4.gt.0.0d0) then
220 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
221 & 'between contact pairs of peptide groups'
222 write (iout,'(2(a,f5.3/))')
223 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
224 & 'Range of quenching the correlation terms:',2*delt_corr
225 else if (wcorr.gt.0.0d0) then
226 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
227 & 'between contact pairs of peptide groups'
229 write (iout,'(a,f8.3)')
230 & 'Scaling factor of 1,4 SC-p interactions:',scal14
231 write (iout,'(a,f8.3)')
232 & 'General scaling factor of SC-p interactions:',scalscp
233 r0_corr=cutoff_corr-delt_corr
235 aad(i,1)=scalscp*aad(i,1)
236 aad(i,2)=scalscp*aad(i,2)
237 bad(i,1)=scalscp*bad(i,1)
238 bad(i,2)=scalscp*bad(i,2)
242 print *,'indpdb=',indpdb,' pdbref=',pdbref
244 C Read sequence if not taken from the pdb file.
245 if (iscode.gt.0) then
246 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
248 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
250 C Convert sequence to numeric code
252 itype(i)=rescode(i,sequence(i),iscode)
255 print '(20i4)',(itype(i),i=1,nres)
259 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
261 if (itype(i).eq.ntyp1) then
265 else if (iabs(itype(i+1)).ne.20) then
267 else if (iabs(itype(i)).ne.20) then
274 write (iout,*) "ITEL"
276 write (iout,*) i,itype(i),itel(i)
279 print *,'Call Read_Bridge.'
283 print *,'NNT=',NNT,' NCT=',NCT
284 if (itype(1).eq.ntyp1) nnt=2
285 if (itype(nres).eq.ntyp1) nct=nct-1
286 if (nstart.lt.nnt) nstart=nnt
287 if (nend.gt.nct .or. nend.eq.0) nend=nct
288 write (iout,*) "nstart",nstart," nend",nend
291 c read(inp,'(a)') pdbfile
292 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
293 c open(ipdbin,file=pdbfile,status='old',err=33)
295 c 33 write (iout,'(a)') 'Error opening PDB file.'
298 c print *,'Begin reading pdb data'
300 c print *,'Finished reading pdb data'
301 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
303 c itype_pdb(i)=itype(i)
306 c write (iout,'(a,i3)') 'nsup=',nsup
308 c if (nsup.le.(nct-nnt+1)) then
309 c do i=0,nct-nnt+1-nsup
310 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
316 c & 'Error - sequences to be superposed do not match.'
319 c do i=0,nsup-(nct-nnt+1)
320 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
322 c nstart_sup=nstart_sup+i
328 c & 'Error - sequences to be superposed do not match.'
331 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
332 c & ' nstart_seq=',nstart_seq
336 write (iout,*) "molread: REFSTR",refstr
338 if (.not.pdbref) then
339 call read_angles(inp,*38)
341 38 write (iout,'(a)') 'Error reading reference structure.'
343 call mp_stopall(Error_Msg)
345 stop 'Error reading reference structure'
358 call contact(.true.,ncont_ref,icont_ref)
361 C write (iout,'(/a,i3,a)')
362 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
363 write (iout,'(20i4)') (iss(i),i=1,ns)
365 write(iout,*)"Running with dynamic disulfide-bond formation"
367 write (iout,'(/a/)') 'Pre-formed links are:'
373 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
374 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
380 if (ns.gt.0.and.dyn_ss) then
384 forcon(i-nss)=forcon(i)
391 dyn_ss_mask(iss(i))=.true.
396 c-----------------------------------------------------------------------------
397 logical function seq_comp(itypea,itypeb,length)
399 integer length,itypea(length),itypeb(length)
402 if (itypea(i).ne.itypeb(i)) then
410 c-----------------------------------------------------------------------------
411 subroutine read_bridge
412 C Read information about disulfide bridges.
415 include 'COMMON.IOUNITS'
418 include 'COMMON.INTERACT'
419 include 'COMMON.LOCAL'
420 include 'COMMON.NAMES'
421 include 'COMMON.CHAIN'
422 include 'COMMON.FFIELD'
423 include 'COMMON.SBRIDGE'
424 include 'COMMON.HEADER'
425 include 'COMMON.CONTROL'
426 include 'COMMON.TIME1'
428 include 'COMMON.INFO'
431 C Read bridging residues.
432 read (inp,*) ns,(iss(i),i=1,ns)
434 C Check whether the specified bridging residues are cystines.
436 if (itype(iss(i)).ne.1) then
437 write (iout,'(2a,i3,a)')
438 & 'Do you REALLY think that the residue ',
439 & restyp(itype(iss(i))),i,
440 & ' can form a disulfide bridge?!!!'
441 write (*,'(2a,i3,a)')
442 & 'Do you REALLY think that the residue ',
443 & restyp(itype(iss(i))),i,
444 & ' can form a disulfide bridge?!!!'
446 call mp_stopall(error_msg)
452 C Read preformed bridges.
454 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
457 C Check if the residues involved in bridges are in the specified list of
461 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
462 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
463 write (iout,'(a,i3,a)') 'Disulfide pair',i,
464 & ' contains residues present in other pairs.'
465 write (*,'(a,i3,a)') 'Disulfide pair',i,
466 & ' contains residues present in other pairs.'
468 call mp_stopall(error_msg)
475 if (ihpb(i).eq.iss(j)) goto 10
477 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
480 if (jhpb(i).eq.iss(j)) goto 20
482 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
495 c----------------------------------------------------------------------------
496 subroutine read_angles(kanal,*)
501 include 'COMMON.CHAIN'
502 include 'COMMON.IOUNITS'
504 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
505 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
506 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
507 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
509 theta(i)=deg2rad*theta(i)
510 phi(i)=deg2rad*phi(i)
511 alph(i)=deg2rad*alph(i)
512 omeg(i)=deg2rad*omeg(i)
517 c----------------------------------------------------------------------------
518 subroutine reada(rekord,lancuch,wartosc,default)
520 character*(*) rekord,lancuch
521 double precision wartosc,default
524 iread=index(rekord,lancuch)
529 iread=iread+ilen(lancuch)+1
530 read (rekord(iread:),*) wartosc
533 c----------------------------------------------------------------------------
534 subroutine multreada(rekord,lancuch,tablica,dim,default)
537 double precision tablica(dim),default
538 character*(*) rekord,lancuch
544 iread=index(rekord,lancuch)
545 if (iread.eq.0) return
546 iread=iread+ilen(lancuch)+1
547 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
550 c----------------------------------------------------------------------------
551 subroutine readi(rekord,lancuch,wartosc,default)
553 character*(*) rekord,lancuch
554 integer wartosc,default
557 iread=index(rekord,lancuch)
562 iread=iread+ilen(lancuch)+1
563 read (rekord(iread:),*) wartosc
566 c----------------------------------------------------------------------------
567 subroutine card_concat(card)
569 include 'COMMON.IOUNITS'
571 character*80 karta,ucase
573 read (inp,'(a)') karta
576 do while (karta(80:80).eq.'&')
577 card=card(:ilen(card)+1)//karta(:79)
578 read (inp,'(a)') karta
581 card=card(:ilen(card)+1)//karta
584 c----------------------------------------------------------------------------
593 include 'COMMON.IOUNITS'
594 include 'COMMON.CONTROL'
595 integer lenpre,lenpot,ilen
597 character*16 cformat,cprint
599 integer lenint,lenout
600 call getenv('INPUT',prefix)
601 call getenv('OUTPUT',prefout)
602 call getenv('INTIN',prefintin)
603 call getenv('COORD',cformat)
604 call getenv('PRINTCOOR',cprint)
605 call getenv('SCRATCHDIR',scratchdir)
608 if (index(ucase(cformat),'CX').gt.0) then
615 lenint=ilen(prefintin)
616 C Get the names and open the input files
617 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
619 write (liczba,'(bz,i3.3)') me
620 outname=prefout(:lenout)//'_clust.out_'//liczba
622 outname=prefout(:lenout)//'_clust.out'
625 intinname=prefintin(:lenint)//'.bx'
626 else if (from_cx) then
627 intinname=prefintin(:lenint)//'.cx'
629 intinname=prefintin(:lenint)//'.int'
631 rmsname=prefintin(:lenint)//'.rms'
632 open (jplot,file=prefout(:ilen(prefout))//'.tex',
634 open (jrms,file=rmsname,status='unknown')
635 open(iout,file=outname,status='unknown')
636 C Get parameter filenames and open the parameter files.
637 call getenv('BONDPAR',bondname)
638 open (ibond,file=bondname,status='old')
639 call getenv('THETPAR',thetname)
640 open (ithep,file=thetname,status='old')
641 call getenv('ROTPAR',rotname)
642 open (irotam,file=rotname,status='old')
643 call getenv('TORPAR',torname)
644 open (itorp,file=torname,status='old')
645 call getenv('TORDPAR',tordname)
646 open (itordp,file=tordname,status='old')
647 call getenv('FOURIER',fouriername)
648 open (ifourier,file=fouriername,status='old')
649 call getenv('ELEPAR',elename)
650 open (ielep,file=elename,status='old')
651 call getenv('SIDEPAR',sidename)
652 open (isidep,file=sidename,status='old')
653 call getenv('SIDEP',sidepname)
654 open (isidep1,file=sidepname,status="old")
655 call getenv('SCCORPAR',sccorname)
656 open (isccor,file=sccorname,status="old")
659 C 8/9/01 In the newest version SCp interaction constants are read from a file
660 C Use -DOLDSCP to use hard-coded constants instead.
662 call getenv('SCPPAR',scpname)
663 open (iscpp,file=scpname,status='old')