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 character*320 controlcard,ucase
23 read (INP,'(a80)') titel
24 call card_concat(controlcard)
26 call readi(controlcard,'NRES',nres,0)
27 call readi(controlcard,'RESCALE',rescale_mode,2)
28 call readi(controlcard,'PDBOUT',outpdb,0)
29 call readi(controlcard,'MOL2OUT',outmol2,0)
30 refstr=(index(controlcard,'REFSTR').gt.0)
31 write (iout,*) "REFSTR",refstr
32 pdbref=(index(controlcard,'PDBREF').gt.0)
33 iscode=index(controlcard,'ONE_LETTER')
34 tree=(index(controlcard,'MAKE_TREE').gt.0)
35 min_var=(index(controlcard,'MINVAR').gt.0)
36 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
37 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
38 call readi(controlcard,'NCUT',ncut,1)
39 call readi(controlcard,'NSTART',nstart,0)
40 call readi(controlcard,'NEND',nend,0)
41 call reada(controlcard,'ECUT',ecut,10.0d0)
42 call reada(controlcard,'PROB',prob_limit,0.99d0)
43 write (iout,*) "Probability limit",prob_limit
44 lgrp=(index(controlcard,'LGRP').gt.0)
45 caonly=(index(controlcard,'CA_ONLY').gt.0)
46 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
47 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
48 call readi(controlcard,'IOPT',iopt,2)
49 lside = index(controlcard,"SIDE").gt.0
50 efree = index(controlcard,"EFREE").gt.0
51 call readi(controlcard,'NTEMP',nT,1)
52 write (iout,*) "nT",nT
53 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
54 write (iout,*) "nT",nT
55 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
57 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
59 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
60 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
61 lprint_int=index(controlcard,"PRINT_INT") .gt.0
65 c--------------------------------------------------------------------------
68 C Read molecular data.
72 include 'COMMON.IOUNITS'
75 include 'COMMON.INTERACT'
76 include 'COMMON.LOCAL'
77 include 'COMMON.NAMES'
78 include 'COMMON.CHAIN'
79 include 'COMMON.FFIELD'
80 include 'COMMON.SBRIDGE'
81 include 'COMMON.HEADER'
82 include 'COMMON.CONTROL'
83 include 'COMMON.CONTACTS'
84 include 'COMMON.TIME1'
88 character*4 sequence(maxres)
89 character*800 weightcard
91 double precision x(maxvar)
92 integer itype_pdb(maxres)
98 C Read weights of the subsequent energy terms.
99 call card_concat(weightcard)
100 call reada(weightcard,'WSC',wsc,1.0d0)
101 call reada(weightcard,'WLONG',wsc,wsc)
102 call reada(weightcard,'WSCP',wscp,1.0d0)
103 call reada(weightcard,'WELEC',welec,1.0D0)
104 call reada(weightcard,'WVDWPP',wvdwpp,welec)
105 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
106 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
107 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
108 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
109 call reada(weightcard,'WTURN3',wturn3,1.0D0)
110 call reada(weightcard,'WTURN4',wturn4,1.0D0)
111 call reada(weightcard,'WTURN6',wturn6,1.0D0)
112 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
113 call reada(weightcard,'WBOND',wbond,1.0D0)
114 call reada(weightcard,'WTOR',wtor,1.0D0)
115 call reada(weightcard,'WTORD',wtor_d,1.0D0)
116 call reada(weightcard,'WANG',wang,1.0D0)
117 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
118 call reada(weightcard,'SCAL14',scal14,0.4D0)
119 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
120 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
121 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
122 if (index(weightcard,'SOFT').gt.0) ipot=6
123 C 12/1/95 Added weight for the multi-body term WCORR
124 call reada(weightcard,'WCORRH',wcorr,1.0D0)
125 if (wcorr4.gt.0.0d0) wcorr=wcorr4
144 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
145 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
147 10 format (/'Energy-term weights (unscaled):'//
148 & 'WSCC= ',f10.6,' (SC-SC)'/
149 & 'WSCP= ',f10.6,' (SC-p)'/
150 & 'WELEC= ',f10.6,' (p-p electr)'/
151 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
152 & 'WBOND= ',f10.6,' (stretching)'/
153 & 'WANG= ',f10.6,' (bending)'/
154 & 'WSCLOC= ',f10.6,' (SC local)'/
155 & 'WTOR= ',f10.6,' (torsional)'/
156 & 'WTORD= ',f10.6,' (double torsional)'/
157 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
158 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
159 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
160 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
161 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
162 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
163 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
164 & 'WTURN6= ',f10.6,' (turns, 6th order)')
165 if (wcorr4.gt.0.0d0) then
166 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
167 & 'between contact pairs of peptide groups'
168 write (iout,'(2(a,f5.3/))')
169 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
170 & 'Range of quenching the correlation terms:',2*delt_corr
171 else if (wcorr.gt.0.0d0) then
172 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
173 & 'between contact pairs of peptide groups'
175 write (iout,'(a,f8.3)')
176 & 'Scaling factor of 1,4 SC-p interactions:',scal14
177 write (iout,'(a,f8.3)')
178 & 'General scaling factor of SC-p interactions:',scalscp
179 r0_corr=cutoff_corr-delt_corr
181 aad(i,1)=scalscp*aad(i,1)
182 aad(i,2)=scalscp*aad(i,2)
183 bad(i,1)=scalscp*bad(i,1)
184 bad(i,2)=scalscp*bad(i,2)
188 print *,'indpdb=',indpdb,' pdbref=',pdbref
190 C Read sequence if not taken from the pdb file.
191 if (iscode.gt.0) then
192 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
194 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
196 C Convert sequence to numeric code
198 itype(i)=rescode(i,sequence(i),iscode)
201 print '(20i4)',(itype(i),i=1,nres)
205 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
207 if (itype(i).eq.21) then
211 else if (itype(i+1).ne.20) then
213 else if (itype(i).ne.20) then
220 write (iout,*) "ITEL"
222 write (iout,*) i,itype(i),itel(i)
225 print *,'Call Read_Bridge.'
229 print *,'NNT=',NNT,' NCT=',NCT
230 if (itype(1).eq.21) nnt=2
231 if (itype(nres).eq.21) nct=nct-1
232 if (nstart.lt.nnt) nstart=nnt
233 if (nend.gt.nct .or. nend.eq.0) nend=nct
234 write (iout,*) "nstart",nstart," nend",nend
237 c read(inp,'(a)') pdbfile
238 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
239 c open(ipdbin,file=pdbfile,status='old',err=33)
241 c 33 write (iout,'(a)') 'Error opening PDB file.'
244 c print *,'Begin reading pdb data'
246 c print *,'Finished reading pdb data'
247 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
249 c itype_pdb(i)=itype(i)
252 c write (iout,'(a,i3)') 'nsup=',nsup
254 c if (nsup.le.(nct-nnt+1)) then
255 c do i=0,nct-nnt+1-nsup
256 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
262 c & 'Error - sequences to be superposed do not match.'
265 c do i=0,nsup-(nct-nnt+1)
266 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
268 c nstart_sup=nstart_sup+i
274 c & 'Error - sequences to be superposed do not match.'
277 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
278 c & ' nstart_seq=',nstart_seq
282 write (iout,*) "molread: REFSTR",refstr
284 if (.not.pdbref) then
285 call read_angles(inp,*38)
287 38 write (iout,'(a)') 'Error reading reference structure.'
289 call mp_stopall(Error_Msg)
291 stop 'Error reading reference structure'
303 call contact(.true.,ncont_ref,icont_ref)
307 c-----------------------------------------------------------------------------
308 logical function seq_comp(itypea,itypeb,length)
310 integer length,itypea(length),itypeb(length)
313 if (itypea(i).ne.itypeb(i)) then
321 c-----------------------------------------------------------------------------
322 subroutine read_bridge
323 C Read information about disulfide bridges.
326 include 'COMMON.IOUNITS'
329 include 'COMMON.INTERACT'
330 include 'COMMON.LOCAL'
331 include 'COMMON.NAMES'
332 include 'COMMON.CHAIN'
333 include 'COMMON.FFIELD'
334 include 'COMMON.SBRIDGE'
335 include 'COMMON.HEADER'
336 include 'COMMON.CONTROL'
337 include 'COMMON.TIME1'
339 include 'COMMON.INFO'
342 C Read bridging residues.
343 read (inp,*) ns,(iss(i),i=1,ns)
345 C Check whether the specified bridging residues are cystines.
347 if (itype(iss(i)).ne.1) then
348 write (iout,'(2a,i3,a)')
349 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
350 & ' can form a disulfide bridge?!!!'
351 write (*,'(2a,i3,a)')
352 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
353 & ' can form a disulfide bridge?!!!'
355 call mp_stopall(error_msg)
361 C Read preformed bridges.
363 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
366 C Check if the residues involved in bridges are in the specified list of
370 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
371 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
372 write (iout,'(a,i3,a)') 'Disulfide pair',i,
373 & ' contains residues present in other pairs.'
374 write (*,'(a,i3,a)') 'Disulfide pair',i,
375 & ' contains residues present in other pairs.'
377 call mp_stopall(error_msg)
384 if (ihpb(i).eq.iss(j)) goto 10
386 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
389 if (jhpb(i).eq.iss(j)) goto 20
391 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
404 c----------------------------------------------------------------------------
405 subroutine read_angles(kanal,*)
410 include 'COMMON.CHAIN'
411 include 'COMMON.IOUNITS'
413 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
414 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
415 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
416 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
418 theta(i)=deg2rad*theta(i)
419 phi(i)=deg2rad*phi(i)
420 alph(i)=deg2rad*alph(i)
421 omeg(i)=deg2rad*omeg(i)
426 c----------------------------------------------------------------------------
427 subroutine reada(rekord,lancuch,wartosc,default)
429 character*(*) rekord,lancuch
430 double precision wartosc,default
433 iread=index(rekord,lancuch)
438 iread=iread+ilen(lancuch)+1
439 read (rekord(iread:),*) wartosc
442 c----------------------------------------------------------------------------
443 subroutine multreada(rekord,lancuch,tablica,dim,default)
446 double precision tablica(dim),default
447 character*(*) rekord,lancuch
453 iread=index(rekord,lancuch)
454 if (iread.eq.0) return
455 iread=iread+ilen(lancuch)+1
456 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
459 c----------------------------------------------------------------------------
460 subroutine readi(rekord,lancuch,wartosc,default)
462 character*(*) rekord,lancuch
463 integer wartosc,default
466 iread=index(rekord,lancuch)
471 iread=iread+ilen(lancuch)+1
472 read (rekord(iread:),*) wartosc
475 c----------------------------------------------------------------------------
476 subroutine card_concat(card)
478 include 'COMMON.IOUNITS'
480 character*80 karta,ucase
482 read (inp,'(a)') karta
485 do while (karta(80:80).eq.'&')
486 card=card(:ilen(card)+1)//karta(:79)
487 read (inp,'(a)') karta
490 card=card(:ilen(card)+1)//karta
493 c----------------------------------------------------------------------------
502 include 'COMMON.IOUNITS'
503 include 'COMMON.CONTROL'
504 integer lenpre,lenpot,ilen
506 character*16 cformat,cprint
508 integer lenint,lenout
509 call getenv('INPUT',prefix)
510 call getenv('OUTPUT',prefout)
511 call getenv('INTIN',prefintin)
512 call getenv('COORD',cformat)
513 call getenv('PRINTCOOR',cprint)
514 call getenv('SCRATCHDIR',scratchdir)
517 if (index(ucase(cformat),'CX').gt.0) then
524 lenint=ilen(prefintin)
525 C Get the names and open the input files
526 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
528 write (liczba,'(bz,i3.3)') me
529 outname=prefout(:lenout)//'_clust.out_'//liczba
531 outname=prefout(:lenout)//'_clust.out'
534 intinname=prefintin(:lenint)//'.bx'
535 else if (from_cx) then
536 intinname=prefintin(:lenint)//'.cx'
538 intinname=prefintin(:lenint)//'.int'
540 rmsname=prefintin(:lenint)//'.rms'
541 open (jplot,file=prefout(:ilen(prefout))//'.tex',
543 open (jrms,file=rmsname,status='unknown')
544 open(iout,file=outname,status='unknown')
545 C Get parameter filenames and open the parameter files.
546 call getenv('BONDPAR',bondname)
547 open (ibond,file=bondname,status='old')
548 call getenv('THETPAR',thetname)
549 open (ithep,file=thetname,status='old')
550 call getenv('ROTPAR',rotname)
551 open (irotam,file=rotname,status='old')
552 call getenv('TORPAR',torname)
553 open (itorp,file=torname,status='old')
554 call getenv('TORDPAR',tordname)
555 open (itordp,file=tordname,status='old')
556 call getenv('FOURIER',fouriername)
557 open (ifourier,file=fouriername,status='old')
558 call getenv('ELEPAR',elename)
559 open (ielep,file=elename,status='old')
560 call getenv('SIDEPAR',sidename)
561 open (isidep,file=sidename,status='old')
562 call getenv('SIDEP',sidepname)
563 open (isidep1,file=sidepname,status="old")
564 call getenv('SCCORPAR',sccorname)
565 open (isccor,file=sccorname,status="old")
568 C 8/9/01 In the newest version SCp interaction constants are read from a file
569 C Use -DOLDSCP to use hard-coded constants instead.
571 call getenv('SCPPAR',scpname)
572 open (iscpp,file=scpname,status='old')