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,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)
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 C 12/1/95 Added weight for the multi-body term WCORR
130 call reada(weightcard,'WCORRH',wcorr,1.0D0)
131 if (wcorr4.gt.0.0d0) wcorr=wcorr4
151 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
152 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
153 & wturn4,wturn6,wsccor
154 10 format (/'Energy-term weights (unscaled):'//
155 & 'WSCC= ',f10.6,' (SC-SC)'/
156 & 'WSCP= ',f10.6,' (SC-p)'/
157 & 'WELEC= ',f10.6,' (p-p electr)'/
158 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
159 & 'WBOND= ',f10.6,' (stretching)'/
160 & 'WANG= ',f10.6,' (bending)'/
161 & 'WSCLOC= ',f10.6,' (SC local)'/
162 & 'WTOR= ',f10.6,' (torsional)'/
163 & 'WTORD= ',f10.6,' (double torsional)'/
164 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
165 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
166 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
167 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
168 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
169 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
170 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
171 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
172 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
174 if (wcorr4.gt.0.0d0) then
175 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
176 & 'between contact pairs of peptide groups'
177 write (iout,'(2(a,f5.3/))')
178 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
179 & 'Range of quenching the correlation terms:',2*delt_corr
180 else if (wcorr.gt.0.0d0) then
181 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
182 & 'between contact pairs of peptide groups'
184 write (iout,'(a,f8.3)')
185 & 'Scaling factor of 1,4 SC-p interactions:',scal14
186 write (iout,'(a,f8.3)')
187 & 'General scaling factor of SC-p interactions:',scalscp
188 r0_corr=cutoff_corr-delt_corr
190 aad(i,1)=scalscp*aad(i,1)
191 aad(i,2)=scalscp*aad(i,2)
192 bad(i,1)=scalscp*bad(i,1)
193 bad(i,2)=scalscp*bad(i,2)
197 print *,'indpdb=',indpdb,' pdbref=',pdbref
199 C Read sequence if not taken from the pdb file.
200 if (iscode.gt.0) then
201 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
203 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
205 C Convert sequence to numeric code
207 itype(i)=rescode(i,sequence(i),iscode)
210 print '(20i4)',(itype(i),i=1,nres)
214 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
216 if (itype(i).eq.ntyp1) then
220 else if (iabs(itype(i+1)).ne.20) then
222 else if (iabs(itype(i)).ne.20) then
229 write (iout,*) "ITEL"
231 write (iout,*) i,itype(i),itel(i)
234 print *,'Call Read_Bridge.'
238 print *,'NNT=',NNT,' NCT=',NCT
239 if (itype(1).eq.ntyp1) nnt=2
240 if (itype(nres).eq.ntyp1) nct=nct-1
241 if (nstart.lt.nnt) nstart=nnt
242 if (nend.gt.nct .or. nend.eq.0) nend=nct
243 write (iout,*) "nstart",nstart," nend",nend
246 c read(inp,'(a)') pdbfile
247 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
248 c open(ipdbin,file=pdbfile,status='old',err=33)
250 c 33 write (iout,'(a)') 'Error opening PDB file.'
253 c print *,'Begin reading pdb data'
255 c print *,'Finished reading pdb data'
256 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
258 c itype_pdb(i)=itype(i)
261 c write (iout,'(a,i3)') 'nsup=',nsup
263 c if (nsup.le.(nct-nnt+1)) then
264 c do i=0,nct-nnt+1-nsup
265 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
271 c & 'Error - sequences to be superposed do not match.'
274 c do i=0,nsup-(nct-nnt+1)
275 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
277 c nstart_sup=nstart_sup+i
283 c & 'Error - sequences to be superposed do not match.'
286 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
287 c & ' nstart_seq=',nstart_seq
291 write (iout,*) "molread: REFSTR",refstr
293 if (.not.pdbref) then
294 call read_angles(inp,*38)
296 38 write (iout,'(a)') 'Error reading reference structure.'
298 call mp_stopall(Error_Msg)
300 stop 'Error reading reference structure'
313 call contact(.true.,ncont_ref,icont_ref)
317 c-----------------------------------------------------------------------------
318 logical function seq_comp(itypea,itypeb,length)
320 integer length,itypea(length),itypeb(length)
323 if (itypea(i).ne.itypeb(i)) then
331 c-----------------------------------------------------------------------------
332 subroutine read_bridge
333 C Read information about disulfide bridges.
336 include 'COMMON.IOUNITS'
339 include 'COMMON.INTERACT'
340 include 'COMMON.LOCAL'
341 include 'COMMON.NAMES'
342 include 'COMMON.CHAIN'
343 include 'COMMON.FFIELD'
344 include 'COMMON.SBRIDGE'
345 include 'COMMON.HEADER'
346 include 'COMMON.CONTROL'
347 include 'COMMON.TIME1'
349 include 'COMMON.INFO'
352 C Read bridging residues.
353 read (inp,*) ns,(iss(i),i=1,ns)
355 C Check whether the specified bridging residues are cystines.
357 if (itype(iss(i)).ne.1) then
358 write (iout,'(2a,i3,a)')
359 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
360 & ' can form a disulfide bridge?!!!'
361 write (*,'(2a,i3,a)')
362 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
363 & ' can form a disulfide bridge?!!!'
365 call mp_stopall(error_msg)
371 C Read preformed bridges.
373 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
376 C Check if the residues involved in bridges are in the specified list of
380 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
381 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
382 write (iout,'(a,i3,a)') 'Disulfide pair',i,
383 & ' contains residues present in other pairs.'
384 write (*,'(a,i3,a)') 'Disulfide pair',i,
385 & ' contains residues present in other pairs.'
387 call mp_stopall(error_msg)
394 if (ihpb(i).eq.iss(j)) goto 10
396 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
399 if (jhpb(i).eq.iss(j)) goto 20
401 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
414 c----------------------------------------------------------------------------
415 subroutine read_angles(kanal,*)
420 include 'COMMON.CHAIN'
421 include 'COMMON.IOUNITS'
423 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
424 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
425 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
426 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
428 theta(i)=deg2rad*theta(i)
429 phi(i)=deg2rad*phi(i)
430 alph(i)=deg2rad*alph(i)
431 omeg(i)=deg2rad*omeg(i)
436 c----------------------------------------------------------------------------
437 subroutine reada(rekord,lancuch,wartosc,default)
439 character*(*) rekord,lancuch
440 double precision wartosc,default
443 iread=index(rekord,lancuch)
448 iread=iread+ilen(lancuch)+1
449 read (rekord(iread:),*) wartosc
452 c----------------------------------------------------------------------------
453 subroutine multreada(rekord,lancuch,tablica,dim,default)
456 double precision tablica(dim),default
457 character*(*) rekord,lancuch
463 iread=index(rekord,lancuch)
464 if (iread.eq.0) return
465 iread=iread+ilen(lancuch)+1
466 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
469 c----------------------------------------------------------------------------
470 subroutine readi(rekord,lancuch,wartosc,default)
472 character*(*) rekord,lancuch
473 integer wartosc,default
476 iread=index(rekord,lancuch)
481 iread=iread+ilen(lancuch)+1
482 read (rekord(iread:),*) wartosc
485 c----------------------------------------------------------------------------
486 subroutine card_concat(card)
488 include 'COMMON.IOUNITS'
490 character*80 karta,ucase
492 read (inp,'(a)') karta
495 do while (karta(80:80).eq.'&')
496 card=card(:ilen(card)+1)//karta(:79)
497 read (inp,'(a)') karta
500 card=card(:ilen(card)+1)//karta
503 c----------------------------------------------------------------------------
512 include 'COMMON.IOUNITS'
513 include 'COMMON.CONTROL'
514 integer lenpre,lenpot,ilen
516 character*16 cformat,cprint
518 integer lenint,lenout
519 call getenv('INPUT',prefix)
520 call getenv('OUTPUT',prefout)
521 call getenv('INTIN',prefintin)
522 call getenv('COORD',cformat)
523 call getenv('PRINTCOOR',cprint)
524 call getenv('SCRATCHDIR',scratchdir)
527 if (index(ucase(cformat),'CX').gt.0) then
534 lenint=ilen(prefintin)
535 C Get the names and open the input files
536 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
538 write (liczba,'(bz,i3.3)') me
539 outname=prefout(:lenout)//'_clust.out_'//liczba
541 outname=prefout(:lenout)//'_clust.out'
544 intinname=prefintin(:lenint)//'.bx'
545 else if (from_cx) then
546 intinname=prefintin(:lenint)//'.cx'
548 intinname=prefintin(:lenint)//'.int'
550 rmsname=prefintin(:lenint)//'.rms'
551 open (jplot,file=prefout(:ilen(prefout))//'.tex',
553 open (jrms,file=rmsname,status='unknown')
554 open(iout,file=outname,status='unknown')
555 C Get parameter filenames and open the parameter files.
556 call getenv('BONDPAR',bondname)
557 open (ibond,file=bondname,status='old')
558 call getenv('THETPAR',thetname)
559 open (ithep,file=thetname,status='old')
560 call getenv('ROTPAR',rotname)
561 open (irotam,file=rotname,status='old')
562 call getenv('TORPAR',torname)
563 open (itorp,file=torname,status='old')
564 call getenv('TORDPAR',tordname)
565 open (itordp,file=tordname,status='old')
566 call getenv('FOURIER',fouriername)
567 open (ifourier,file=fouriername,status='old')
568 call getenv('ELEPAR',elename)
569 open (ielep,file=elename,status='old')
570 call getenv('SIDEPAR',sidename)
571 open (isidep,file=sidename,status='old')
572 call getenv('SIDEP',sidepname)
573 open (isidep1,file=sidepname,status="old")
574 call getenv('SCCORPAR',sccorname)
575 open (isccor,file=sccorname,status="old")
578 C 8/9/01 In the newest version SCp interaction constants are read from a file
579 C Use -DOLDSCP to use hard-coded constants instead.
581 call getenv('SCPPAR',scpname)
582 open (iscpp,file=scpname,status='old')