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,'WBOND',wbond,1.0D0)
119 call reada(weightcard,'WTOR',wtor,1.0D0)
120 call reada(weightcard,'WTORD',wtor_d,1.0D0)
121 call reada(weightcard,'WANG',wang,1.0D0)
122 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
123 call reada(weightcard,'SCAL14',scal14,0.4D0)
124 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
125 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
126 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
127 if (index(weightcard,'SOFT').gt.0) ipot=6
128 C 12/1/95 Added weight for the multi-body term WCORR
129 call reada(weightcard,'WCORRH',wcorr,1.0D0)
130 if (wcorr4.gt.0.0d0) wcorr=wcorr4
149 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
150 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
152 10 format (/'Energy-term weights (unscaled):'//
153 & 'WSCC= ',f10.6,' (SC-SC)'/
154 & 'WSCP= ',f10.6,' (SC-p)'/
155 & 'WELEC= ',f10.6,' (p-p electr)'/
156 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
157 & 'WBOND= ',f10.6,' (stretching)'/
158 & 'WANG= ',f10.6,' (bending)'/
159 & 'WSCLOC= ',f10.6,' (SC local)'/
160 & 'WTOR= ',f10.6,' (torsional)'/
161 & 'WTORD= ',f10.6,' (double torsional)'/
162 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
163 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
164 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
165 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
166 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
167 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
168 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
169 & 'WTURN6= ',f10.6,' (turns, 6th order)')
170 if (wcorr4.gt.0.0d0) then
171 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
172 & 'between contact pairs of peptide groups'
173 write (iout,'(2(a,f5.3/))')
174 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
175 & 'Range of quenching the correlation terms:',2*delt_corr
176 else if (wcorr.gt.0.0d0) then
177 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
178 & 'between contact pairs of peptide groups'
180 write (iout,'(a,f8.3)')
181 & 'Scaling factor of 1,4 SC-p interactions:',scal14
182 write (iout,'(a,f8.3)')
183 & 'General scaling factor of SC-p interactions:',scalscp
184 r0_corr=cutoff_corr-delt_corr
186 aad(i,1)=scalscp*aad(i,1)
187 aad(i,2)=scalscp*aad(i,2)
188 bad(i,1)=scalscp*bad(i,1)
189 bad(i,2)=scalscp*bad(i,2)
193 print *,'indpdb=',indpdb,' pdbref=',pdbref
195 C Read sequence if not taken from the pdb file.
196 if (iscode.gt.0) then
197 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
199 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
201 C Convert sequence to numeric code
203 itype(i)=rescode(i,sequence(i),iscode)
206 print '(20i4)',(itype(i),i=1,nres)
210 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
212 if (itype(i).eq.21) then
216 else if (itype(i+1).ne.20) then
218 else if (itype(i).ne.20) then
225 write (iout,*) "ITEL"
227 write (iout,*) i,itype(i),itel(i)
230 print *,'Call Read_Bridge.'
234 print *,'NNT=',NNT,' NCT=',NCT
235 if (itype(1).eq.21) nnt=2
236 if (itype(nres).eq.21) nct=nct-1
237 if (nstart.lt.nnt) nstart=nnt
238 if (nend.gt.nct .or. nend.eq.0) nend=nct
239 write (iout,*) "nstart",nstart," nend",nend
242 c read(inp,'(a)') pdbfile
243 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
244 c open(ipdbin,file=pdbfile,status='old',err=33)
246 c 33 write (iout,'(a)') 'Error opening PDB file.'
249 c print *,'Begin reading pdb data'
251 c print *,'Finished reading pdb data'
252 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
254 c itype_pdb(i)=itype(i)
257 c write (iout,'(a,i3)') 'nsup=',nsup
259 c if (nsup.le.(nct-nnt+1)) then
260 c do i=0,nct-nnt+1-nsup
261 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
267 c & 'Error - sequences to be superposed do not match.'
270 c do i=0,nsup-(nct-nnt+1)
271 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
273 c nstart_sup=nstart_sup+i
279 c & 'Error - sequences to be superposed do not match.'
282 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
283 c & ' nstart_seq=',nstart_seq
287 write (iout,*) "molread: REFSTR",refstr
289 if (.not.pdbref) then
290 call read_angles(inp,*38)
292 38 write (iout,'(a)') 'Error reading reference structure.'
294 call mp_stopall(Error_Msg)
296 stop 'Error reading reference structure'
308 call contact(.true.,ncont_ref,icont_ref)
312 c-----------------------------------------------------------------------------
313 logical function seq_comp(itypea,itypeb,length)
315 integer length,itypea(length),itypeb(length)
318 if (itypea(i).ne.itypeb(i)) then
326 c-----------------------------------------------------------------------------
327 subroutine read_bridge
328 C Read information about disulfide bridges.
331 include 'COMMON.IOUNITS'
334 include 'COMMON.INTERACT'
335 include 'COMMON.LOCAL'
336 include 'COMMON.NAMES'
337 include 'COMMON.CHAIN'
338 include 'COMMON.FFIELD'
339 include 'COMMON.SBRIDGE'
340 include 'COMMON.HEADER'
341 include 'COMMON.CONTROL'
342 include 'COMMON.TIME1'
344 include 'COMMON.INFO'
347 C Read bridging residues.
348 read (inp,*) ns,(iss(i),i=1,ns)
350 C Check whether the specified bridging residues are cystines.
352 if (itype(iss(i)).ne.1) then
353 write (iout,'(2a,i3,a)')
354 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
355 & ' can form a disulfide bridge?!!!'
356 write (*,'(2a,i3,a)')
357 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
358 & ' can form a disulfide bridge?!!!'
360 call mp_stopall(error_msg)
366 C Read preformed bridges.
368 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
371 C Check if the residues involved in bridges are in the specified list of
375 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
376 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
377 write (iout,'(a,i3,a)') 'Disulfide pair',i,
378 & ' contains residues present in other pairs.'
379 write (*,'(a,i3,a)') 'Disulfide pair',i,
380 & ' contains residues present in other pairs.'
382 call mp_stopall(error_msg)
389 if (ihpb(i).eq.iss(j)) goto 10
391 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
394 if (jhpb(i).eq.iss(j)) goto 20
396 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
409 c----------------------------------------------------------------------------
410 subroutine read_angles(kanal,*)
415 include 'COMMON.CHAIN'
416 include 'COMMON.IOUNITS'
418 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
419 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
420 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
421 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
423 theta(i)=deg2rad*theta(i)
424 phi(i)=deg2rad*phi(i)
425 alph(i)=deg2rad*alph(i)
426 omeg(i)=deg2rad*omeg(i)
431 c----------------------------------------------------------------------------
432 subroutine reada(rekord,lancuch,wartosc,default)
434 character*(*) rekord,lancuch
435 double precision wartosc,default
438 iread=index(rekord,lancuch)
443 iread=iread+ilen(lancuch)+1
444 read (rekord(iread:),*) wartosc
447 c----------------------------------------------------------------------------
448 subroutine multreada(rekord,lancuch,tablica,dim,default)
451 double precision tablica(dim),default
452 character*(*) rekord,lancuch
458 iread=index(rekord,lancuch)
459 if (iread.eq.0) return
460 iread=iread+ilen(lancuch)+1
461 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
464 c----------------------------------------------------------------------------
465 subroutine readi(rekord,lancuch,wartosc,default)
467 character*(*) rekord,lancuch
468 integer wartosc,default
471 iread=index(rekord,lancuch)
476 iread=iread+ilen(lancuch)+1
477 read (rekord(iread:),*) wartosc
480 c----------------------------------------------------------------------------
481 subroutine card_concat(card)
483 include 'COMMON.IOUNITS'
485 character*80 karta,ucase
487 read (inp,'(a)') karta
490 do while (karta(80:80).eq.'&')
491 card=card(:ilen(card)+1)//karta(:79)
492 read (inp,'(a)') karta
495 card=card(:ilen(card)+1)//karta
498 c----------------------------------------------------------------------------
507 include 'COMMON.IOUNITS'
508 include 'COMMON.CONTROL'
509 integer lenpre,lenpot,ilen
511 character*16 cformat,cprint
513 integer lenint,lenout
514 call getenv('INPUT',prefix)
515 call getenv('OUTPUT',prefout)
516 call getenv('INTIN',prefintin)
517 call getenv('COORD',cformat)
518 call getenv('PRINTCOOR',cprint)
519 call getenv('SCRATCHDIR',scratchdir)
522 if (index(ucase(cformat),'CX').gt.0) then
529 lenint=ilen(prefintin)
530 C Get the names and open the input files
531 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
533 write (liczba,'(bz,i3.3)') me
534 outname=prefout(:lenout)//'_clust.out_'//liczba
536 outname=prefout(:lenout)//'_clust.out'
539 intinname=prefintin(:lenint)//'.bx'
540 else if (from_cx) then
541 intinname=prefintin(:lenint)//'.cx'
543 intinname=prefintin(:lenint)//'.int'
545 rmsname=prefintin(:lenint)//'.rms'
546 open (jplot,file=prefout(:ilen(prefout))//'.tex',
548 open (jrms,file=rmsname,status='unknown')
549 open(iout,file=outname,status='unknown')
550 C Get parameter filenames and open the parameter files.
551 call getenv('BONDPAR',bondname)
552 open (ibond,file=bondname,status='old')
553 call getenv('THETPAR',thetname)
554 open (ithep,file=thetname,status='old')
555 call getenv('ROTPAR',rotname)
556 open (irotam,file=rotname,status='old')
557 call getenv('TORPAR',torname)
558 open (itorp,file=torname,status='old')
559 call getenv('TORDPAR',tordname)
560 open (itordp,file=tordname,status='old')
561 call getenv('FOURIER',fouriername)
562 open (ifourier,file=fouriername,status='old')
563 call getenv('ELEPAR',elename)
564 open (ielep,file=elename,status='old')
565 call getenv('SIDEPAR',sidename)
566 open (isidep,file=sidename,status='old')
567 call getenv('SIDEP',sidepname)
568 open (isidep1,file=sidepname,status="old")
569 call getenv('SCCORPAR',sccorname)
570 open (isccor,file=sccorname,status="old")
573 C 8/9/01 In the newest version SCp interaction constants are read from a file
574 C Use -DOLDSCP to use hard-coded constants instead.
576 call getenv('SCPPAR',scpname)
577 open (iscpp,file=scpname,status='old')