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
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,0)
42 call readi(controlcard,'NCLUST',nclust,5)
43 call readi(controlcard,'SYM',symetr,1)
44 write (iout,*) 'sym', symetr
45 call readi(controlcard,'NSTART',nstart,0)
46 call readi(controlcard,'NEND',nend,0)
47 call reada(controlcard,'ECUT',ecut,10.0d0)
48 call reada(controlcard,'PROB',prob_limit,0.99d0)
49 write (iout,*) "Probability limit",prob_limit
50 lgrp=(index(controlcard,'LGRP').gt.0)
51 caonly=(index(controlcard,'CA_ONLY').gt.0)
52 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
54 & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
55 call readi(controlcard,'IOPT',iopt,2)
56 lside = index(controlcard,"SIDE").gt.0
57 efree = index(controlcard,"EFREE").gt.0
58 call readi(controlcard,'NTEMP',nT,1)
59 write (iout,*) "nT",nT
60 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
61 write (iout,*) "nT",nT
62 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
64 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
66 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
67 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
68 lprint_int=index(controlcard,"PRINT_INT") .gt.0
72 c--------------------------------------------------------------------------
75 C Read molecular data.
79 include 'COMMON.IOUNITS'
82 include 'COMMON.INTERACT'
83 include 'COMMON.LOCAL'
84 include 'COMMON.NAMES'
85 include 'COMMON.CHAIN'
86 include 'COMMON.FFIELD'
87 include 'COMMON.SBRIDGE'
88 include 'COMMON.HEADER'
89 include 'COMMON.CONTROL'
90 include 'COMMON.CONTACTS'
91 include 'COMMON.TIME1'
95 character*4 sequence(maxres)
96 character*800 weightcard
98 double precision x(maxvar)
99 integer itype_pdb(maxres)
105 C Read weights of the subsequent energy terms.
106 call card_concat(weightcard)
107 call reada(weightcard,'WSC',wsc,1.0d0)
108 call reada(weightcard,'WLONG',wsc,wsc)
109 call reada(weightcard,'WSCP',wscp,1.0d0)
110 call reada(weightcard,'WELEC',welec,1.0D0)
111 call reada(weightcard,'WVDWPP',wvdwpp,welec)
112 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
113 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
114 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
115 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
116 call reada(weightcard,'WTURN3',wturn3,1.0D0)
117 call reada(weightcard,'WTURN4',wturn4,1.0D0)
118 call reada(weightcard,'WTURN6',wturn6,1.0D0)
119 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
120 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
121 call reada(weightcard,'WBOND',wbond,1.0D0)
122 call reada(weightcard,'WTOR',wtor,1.0D0)
123 call reada(weightcard,'WTORD',wtor_d,1.0D0)
124 call reada(weightcard,'WANG',wang,1.0D0)
125 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
126 call reada(weightcard,'SCAL14',scal14,0.4D0)
127 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
128 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
129 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
130 if (index(weightcard,'SOFT').gt.0) ipot=6
131 C 12/1/95 Added weight for the multi-body term WCORR
132 call reada(weightcard,'WCORRH',wcorr,1.0D0)
133 if (wcorr4.gt.0.0d0) wcorr=wcorr4
153 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
154 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
155 & wturn4,wturn6,wsccor
156 10 format (/'Energy-term weights (unscaled):'//
157 & 'WSCC= ',f10.6,' (SC-SC)'/
158 & 'WSCP= ',f10.6,' (SC-p)'/
159 & 'WELEC= ',f10.6,' (p-p electr)'/
160 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
161 & 'WBOND= ',f10.6,' (stretching)'/
162 & 'WANG= ',f10.6,' (bending)'/
163 & 'WSCLOC= ',f10.6,' (SC local)'/
164 & 'WTOR= ',f10.6,' (torsional)'/
165 & 'WTORD= ',f10.6,' (double torsional)'/
166 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
167 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
168 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
169 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
170 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
171 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
172 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
173 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
174 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
176 if (wcorr4.gt.0.0d0) then
177 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
178 & 'between contact pairs of peptide groups'
179 write (iout,'(2(a,f5.3/))')
180 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
181 & 'Range of quenching the correlation terms:',2*delt_corr
182 else if (wcorr.gt.0.0d0) then
183 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
184 & 'between contact pairs of peptide groups'
186 write (iout,'(a,f8.3)')
187 & 'Scaling factor of 1,4 SC-p interactions:',scal14
188 write (iout,'(a,f8.3)')
189 & 'General scaling factor of SC-p interactions:',scalscp
190 r0_corr=cutoff_corr-delt_corr
192 aad(i,1)=scalscp*aad(i,1)
193 aad(i,2)=scalscp*aad(i,2)
194 bad(i,1)=scalscp*bad(i,1)
195 bad(i,2)=scalscp*bad(i,2)
199 print *,'indpdb=',indpdb,' pdbref=',pdbref
201 C Read sequence if not taken from the pdb file.
202 if (iscode.gt.0) then
203 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
205 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
207 C Convert sequence to numeric code
209 itype(i)=rescode(i,sequence(i),iscode)
211 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
213 & "Glycine is the first full residue, initial dummy deleted"
219 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
221 & "Glycine is the last full residue, terminal dummy deleted"
225 print '(20i4)',(itype(i),i=1,nres)
229 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
231 if (itype(i).eq.21) then
235 else if (itype(i+1).ne.20) then
237 else if (itype(i).ne.20) then
244 write (iout,*) "ITEL"
246 write (iout,*) i,itype(i),itel(i)
249 print *,'Call Read_Bridge.'
253 print *,'NNT=',NNT,' NCT=',NCT
254 if (itype(1).eq.21) nnt=2
255 if (itype(nres).eq.21) nct=nct-1
256 if (nstart.lt.nnt) nstart=nnt
257 if (nend.gt.nct .or. nend.eq.0) nend=nct
258 write (iout,*) "nstart",nstart," nend",nend
261 c read(inp,'(a)') pdbfile
262 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
263 c open(ipdbin,file=pdbfile,status='old',err=33)
265 c 33 write (iout,'(a)') 'Error opening PDB file.'
268 c print *,'Begin reading pdb data'
270 c print *,'Finished reading pdb data'
271 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
273 c itype_pdb(i)=itype(i)
276 c write (iout,'(a,i3)') 'nsup=',nsup
278 c if (nsup.le.(nct-nnt+1)) then
279 c do i=0,nct-nnt+1-nsup
280 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
286 c & 'Error - sequences to be superposed do not match.'
289 c do i=0,nsup-(nct-nnt+1)
290 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
292 c nstart_sup=nstart_sup+i
298 c & 'Error - sequences to be superposed do not match.'
301 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
302 c & ' nstart_seq=',nstart_seq
306 write (iout,*) "molread: REFSTR",refstr
308 if (.not.pdbref) then
309 call read_angles(inp,*38)
311 38 write (iout,'(a)') 'Error reading reference structure.'
313 call mp_stopall(Error_Msg)
315 stop 'Error reading reference structure'
327 call contact(.true.,ncont_ref,icont_ref)
331 c-----------------------------------------------------------------------------
332 logical function seq_comp(itypea,itypeb,length)
334 integer length,itypea(length),itypeb(length)
337 if (itypea(i).ne.itypeb(i)) then
345 c-----------------------------------------------------------------------------
346 subroutine read_bridge
347 C Read information about disulfide bridges.
350 include 'COMMON.IOUNITS'
353 include 'COMMON.INTERACT'
354 include 'COMMON.LOCAL'
355 include 'COMMON.NAMES'
356 include 'COMMON.CHAIN'
357 include 'COMMON.FFIELD'
358 include 'COMMON.SBRIDGE'
359 include 'COMMON.HEADER'
360 include 'COMMON.CONTROL'
361 include 'COMMON.TIME1'
363 include 'COMMON.INFO'
366 C Read bridging residues.
367 read (inp,*) ns,(iss(i),i=1,ns)
369 C Check whether the specified bridging residues are cystines.
371 if (itype(iss(i)).ne.1) then
372 write (iout,'(2a,i3,a)')
373 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
374 & ' can form a disulfide bridge?!!!'
375 write (*,'(2a,i3,a)')
376 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
377 & ' can form a disulfide bridge?!!!'
379 call mp_stopall(error_msg)
385 C Read preformed bridges.
387 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
390 C Check if the residues involved in bridges are in the specified list of
394 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
395 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
396 write (iout,'(a,i3,a)') 'Disulfide pair',i,
397 & ' contains residues present in other pairs.'
398 write (*,'(a,i3,a)') 'Disulfide pair',i,
399 & ' contains residues present in other pairs.'
401 call mp_stopall(error_msg)
408 if (ihpb(i).eq.iss(j)) goto 10
410 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
413 if (jhpb(i).eq.iss(j)) goto 20
415 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
428 c----------------------------------------------------------------------------
429 subroutine read_angles(kanal,*)
434 include 'COMMON.CHAIN'
435 include 'COMMON.IOUNITS'
437 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
438 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
439 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
440 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
442 theta(i)=deg2rad*theta(i)
443 phi(i)=deg2rad*phi(i)
444 alph(i)=deg2rad*alph(i)
445 omeg(i)=deg2rad*omeg(i)
450 c----------------------------------------------------------------------------
451 subroutine reada(rekord,lancuch,wartosc,default)
453 character*(*) rekord,lancuch
454 double precision wartosc,default
457 iread=index(rekord,lancuch)
462 iread=iread+ilen(lancuch)+1
463 read (rekord(iread:),*) wartosc
466 c----------------------------------------------------------------------------
467 subroutine multreada(rekord,lancuch,tablica,dim,default)
470 double precision tablica(dim),default
471 character*(*) rekord,lancuch
477 iread=index(rekord,lancuch)
478 if (iread.eq.0) return
479 iread=iread+ilen(lancuch)+1
480 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
483 c----------------------------------------------------------------------------
484 subroutine readi(rekord,lancuch,wartosc,default)
486 character*(*) rekord,lancuch
487 integer wartosc,default
490 iread=index(rekord,lancuch)
495 iread=iread+ilen(lancuch)+1
496 read (rekord(iread:),*) wartosc
499 c----------------------------------------------------------------------------
500 subroutine card_concat(card)
502 include 'COMMON.IOUNITS'
504 character*80 karta,ucase
506 read (inp,'(a)') karta
509 do while (karta(80:80).eq.'&')
510 card=card(:ilen(card)+1)//karta(:79)
511 read (inp,'(a)') karta
514 card=card(:ilen(card)+1)//karta
517 c----------------------------------------------------------------------------
526 include 'COMMON.IOUNITS'
527 include 'COMMON.CONTROL'
528 integer lenpre,lenpot,ilen
530 character*16 cformat,cprint
532 integer lenint,lenout
533 call getenv('INPUT',prefix)
534 call getenv('OUTPUT',prefout)
535 call getenv('INTIN',prefintin)
536 call getenv('COORD',cformat)
537 call getenv('PRINTCOOR',cprint)
538 call getenv('SCRATCHDIR',scratchdir)
541 if (index(ucase(cformat),'CX').gt.0) then
548 lenint=ilen(prefintin)
549 C Get the names and open the input files
550 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
552 write (liczba,'(bz,i3.3)') me
553 outname=prefout(:lenout)//'_clust.out_'//liczba
555 outname=prefout(:lenout)//'_clust.out'
558 intinname=prefintin(:lenint)//'.bx'
559 else if (from_cx) then
560 intinname=prefintin(:lenint)//'.cx'
562 intinname=prefintin(:lenint)//'.int'
564 rmsname=prefintin(:lenint)//'.rms'
565 open (jplot,file=prefout(:ilen(prefout))//'.tex',
567 open (jrms,file=rmsname,status='unknown')
568 open(iout,file=outname,status='unknown')
569 C Get parameter filenames and open the parameter files.
570 call getenv('BONDPAR',bondname)
571 open (ibond,file=bondname,status='old')
572 call getenv('THETPAR',thetname)
573 open (ithep,file=thetname,status='old')
574 call getenv('ROTPAR',rotname)
575 open (irotam,file=rotname,status='old')
576 call getenv('TORPAR',torname)
577 open (itorp,file=torname,status='old')
578 call getenv('TORDPAR',tordname)
579 open (itordp,file=tordname,status='old')
580 call getenv('FOURIER',fouriername)
581 open (ifourier,file=fouriername,status='old')
582 call getenv('ELEPAR',elename)
583 open (ielep,file=elename,status='old')
584 call getenv('SIDEPAR',sidename)
585 open (isidep,file=sidename,status='old')
586 call getenv('SIDEP',sidepname)
587 open (isidep1,file=sidepname,status="old")
588 call getenv('SCCORPAR',sccorname)
589 open (isccor,file=sccorname,status="old")
592 C 8/9/01 In the newest version SCp interaction constants are read from a file
593 C Use -DOLDSCP to use hard-coded constants instead.
595 call getenv('SCPPAR',scpname)
596 open (iscpp,file=scpname,status='old')