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 include "COMMON.SPLITELE"
19 character*320 controlcard,ucase
25 read (INP,'(a80)') titel
26 call card_concat(controlcard)
28 call readi(controlcard,'NRES',nres,0)
29 call readi(controlcard,'RESCALE',rescale_mode,2)
30 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
31 write (iout,*) "DISTCHAINMAX",distchainmax
32 C Reading the dimensions of box in x,y,z coordinates
33 call reada(controlcard,'BOXX',boxxsize,100.0d0)
34 call reada(controlcard,'BOXY',boxysize,100.0d0)
35 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
36 c Cutoff range for interactions
37 call reada(controlcard,"R_CUT",r_cut,15.0d0)
38 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
39 call readi(controlcard,'PDBOUT',outpdb,0)
40 call readi(controlcard,'MOL2OUT',outmol2,0)
41 refstr=(index(controlcard,'REFSTR').gt.0)
42 write (iout,*) "REFSTR",refstr
43 pdbref=(index(controlcard,'PDBREF').gt.0)
44 iscode=index(controlcard,'ONE_LETTER')
45 tree=(index(controlcard,'MAKE_TREE').gt.0)
46 min_var=(index(controlcard,'MINVAR').gt.0)
47 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
48 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
49 call readi(controlcard,'NCUT',ncut,1)
50 call readi(controlcard,'SYM',symetr,1)
51 write (iout,*) 'sym', symetr
52 call readi(controlcard,'NSTART',nstart,0)
53 call readi(controlcard,'NEND',nend,0)
54 call reada(controlcard,'ECUT',ecut,10.0d0)
55 call reada(controlcard,'PROB',prob_limit,0.99d0)
56 write (iout,*) "Probability limit",prob_limit
57 lgrp=(index(controlcard,'LGRP').gt.0)
58 caonly=(index(controlcard,'CA_ONLY').gt.0)
59 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
60 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
61 call readi(controlcard,'IOPT',iopt,2)
62 lside = index(controlcard,"SIDE").gt.0
63 efree = index(controlcard,"EFREE").gt.0
64 call readi(controlcard,'NTEMP',nT,1)
65 write (iout,*) "nT",nT
66 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
67 write (iout,*) "nT",nT
68 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
70 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
72 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
73 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
74 lprint_int=index(controlcard,"PRINT_INT") .gt.0
78 c--------------------------------------------------------------------------
81 C Read molecular data.
85 include 'COMMON.IOUNITS'
88 include 'COMMON.INTERACT'
89 include 'COMMON.LOCAL'
90 include 'COMMON.NAMES'
91 include 'COMMON.CHAIN'
92 include 'COMMON.FFIELD'
93 include 'COMMON.SBRIDGE'
94 include 'COMMON.HEADER'
95 include 'COMMON.CONTROL'
96 include 'COMMON.CONTACTS'
97 include 'COMMON.TIME1'
101 character*4 sequence(maxres)
102 character*800 weightcard
104 double precision x(maxvar)
105 integer itype_pdb(maxres)
111 C Read weights of the subsequent energy terms.
112 call card_concat(weightcard)
113 call reada(weightcard,'WSC',wsc,1.0d0)
114 call reada(weightcard,'WLONG',wsc,wsc)
115 call reada(weightcard,'WSCP',wscp,1.0d0)
116 call reada(weightcard,'WELEC',welec,1.0D0)
117 call reada(weightcard,'WVDWPP',wvdwpp,welec)
118 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
119 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
120 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
121 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
122 call reada(weightcard,'WTURN3',wturn3,1.0D0)
123 call reada(weightcard,'WTURN4',wturn4,1.0D0)
124 call reada(weightcard,'WTURN6',wturn6,1.0D0)
125 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
126 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
127 call reada(weightcard,'WBOND',wbond,1.0D0)
128 call reada(weightcard,'WTOR',wtor,1.0D0)
129 call reada(weightcard,'WTORD',wtor_d,1.0D0)
130 call reada(weightcard,'WANG',wang,1.0D0)
131 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
132 call reada(weightcard,'SCAL14',scal14,0.4D0)
133 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
134 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
135 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
136 if (index(weightcard,'SOFT').gt.0) ipot=6
137 C 12/1/95 Added weight for the multi-body term WCORR
138 call reada(weightcard,'WCORRH',wcorr,1.0D0)
139 if (wcorr4.gt.0.0d0) wcorr=wcorr4
159 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
160 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
161 & wturn4,wturn6,wsccor
162 10 format (/'Energy-term weights (unscaled):'//
163 & 'WSCC= ',f10.6,' (SC-SC)'/
164 & 'WSCP= ',f10.6,' (SC-p)'/
165 & 'WELEC= ',f10.6,' (p-p electr)'/
166 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
167 & 'WBOND= ',f10.6,' (stretching)'/
168 & 'WANG= ',f10.6,' (bending)'/
169 & 'WSCLOC= ',f10.6,' (SC local)'/
170 & 'WTOR= ',f10.6,' (torsional)'/
171 & 'WTORD= ',f10.6,' (double torsional)'/
172 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
173 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
174 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
175 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
176 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
177 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
178 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
179 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
180 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
182 if (wcorr4.gt.0.0d0) then
183 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
184 & 'between contact pairs of peptide groups'
185 write (iout,'(2(a,f5.3/))')
186 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
187 & 'Range of quenching the correlation terms:',2*delt_corr
188 else if (wcorr.gt.0.0d0) then
189 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
190 & 'between contact pairs of peptide groups'
192 write (iout,'(a,f8.3)')
193 & 'Scaling factor of 1,4 SC-p interactions:',scal14
194 write (iout,'(a,f8.3)')
195 & 'General scaling factor of SC-p interactions:',scalscp
196 r0_corr=cutoff_corr-delt_corr
198 aad(i,1)=scalscp*aad(i,1)
199 aad(i,2)=scalscp*aad(i,2)
200 bad(i,1)=scalscp*bad(i,1)
201 bad(i,2)=scalscp*bad(i,2)
205 print *,'indpdb=',indpdb,' pdbref=',pdbref
207 C Read sequence if not taken from the pdb file.
208 if (iscode.gt.0) then
209 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
211 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
213 C Convert sequence to numeric code
215 itype(i)=rescode(i,sequence(i),iscode)
218 print '(20i4)',(itype(i),i=1,nres)
222 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
224 if (itype(i).eq.ntyp1) then
228 else if (iabs(itype(i+1)).ne.20) then
230 else if (iabs(itype(i)).ne.20) then
237 write (iout,*) "ITEL"
239 write (iout,*) i,itype(i),itel(i)
242 print *,'Call Read_Bridge.'
246 print *,'NNT=',NNT,' NCT=',NCT
247 if (itype(1).eq.ntyp1) nnt=2
248 if (itype(nres).eq.ntyp1) nct=nct-1
249 if (nstart.lt.nnt) nstart=nnt
250 if (nend.gt.nct .or. nend.eq.0) nend=nct
251 write (iout,*) "nstart",nstart," nend",nend
254 c read(inp,'(a)') pdbfile
255 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
256 c open(ipdbin,file=pdbfile,status='old',err=33)
258 c 33 write (iout,'(a)') 'Error opening PDB file.'
261 c print *,'Begin reading pdb data'
263 c print *,'Finished reading pdb data'
264 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
266 c itype_pdb(i)=itype(i)
269 c write (iout,'(a,i3)') 'nsup=',nsup
271 c if (nsup.le.(nct-nnt+1)) then
272 c do i=0,nct-nnt+1-nsup
273 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
279 c & 'Error - sequences to be superposed do not match.'
282 c do i=0,nsup-(nct-nnt+1)
283 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
285 c nstart_sup=nstart_sup+i
291 c & 'Error - sequences to be superposed do not match.'
294 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
295 c & ' nstart_seq=',nstart_seq
299 write (iout,*) "molread: REFSTR",refstr
301 if (.not.pdbref) then
302 call read_angles(inp,*38)
304 38 write (iout,'(a)') 'Error reading reference structure.'
306 call mp_stopall(Error_Msg)
308 stop 'Error reading reference structure'
321 call contact(.true.,ncont_ref,icont_ref)
325 c-----------------------------------------------------------------------------
326 logical function seq_comp(itypea,itypeb,length)
328 integer length,itypea(length),itypeb(length)
331 if (itypea(i).ne.itypeb(i)) then
339 c-----------------------------------------------------------------------------
340 subroutine read_bridge
341 C Read information about disulfide bridges.
344 include 'COMMON.IOUNITS'
347 include 'COMMON.INTERACT'
348 include 'COMMON.LOCAL'
349 include 'COMMON.NAMES'
350 include 'COMMON.CHAIN'
351 include 'COMMON.FFIELD'
352 include 'COMMON.SBRIDGE'
353 include 'COMMON.HEADER'
354 include 'COMMON.CONTROL'
355 include 'COMMON.TIME1'
357 include 'COMMON.INFO'
360 C Read bridging residues.
361 read (inp,*) ns,(iss(i),i=1,ns)
363 C Check whether the specified bridging residues are cystines.
365 if (itype(iss(i)).ne.1) then
366 write (iout,'(2a,i3,a)')
367 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
368 & ' can form a disulfide bridge?!!!'
369 write (*,'(2a,i3,a)')
370 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
371 & ' can form a disulfide bridge?!!!'
373 call mp_stopall(error_msg)
379 C Read preformed bridges.
381 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
384 C Check if the residues involved in bridges are in the specified list of
388 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
389 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
390 write (iout,'(a,i3,a)') 'Disulfide pair',i,
391 & ' contains residues present in other pairs.'
392 write (*,'(a,i3,a)') 'Disulfide pair',i,
393 & ' contains residues present in other pairs.'
395 call mp_stopall(error_msg)
402 if (ihpb(i).eq.iss(j)) goto 10
404 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
407 if (jhpb(i).eq.iss(j)) goto 20
409 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
422 c----------------------------------------------------------------------------
423 subroutine read_angles(kanal,*)
428 include 'COMMON.CHAIN'
429 include 'COMMON.IOUNITS'
431 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
432 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
433 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
434 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
436 theta(i)=deg2rad*theta(i)
437 phi(i)=deg2rad*phi(i)
438 alph(i)=deg2rad*alph(i)
439 omeg(i)=deg2rad*omeg(i)
444 c----------------------------------------------------------------------------
445 subroutine reada(rekord,lancuch,wartosc,default)
447 character*(*) rekord,lancuch
448 double precision wartosc,default
451 iread=index(rekord,lancuch)
456 iread=iread+ilen(lancuch)+1
457 read (rekord(iread:),*) wartosc
460 c----------------------------------------------------------------------------
461 subroutine multreada(rekord,lancuch,tablica,dim,default)
464 double precision tablica(dim),default
465 character*(*) rekord,lancuch
471 iread=index(rekord,lancuch)
472 if (iread.eq.0) return
473 iread=iread+ilen(lancuch)+1
474 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
477 c----------------------------------------------------------------------------
478 subroutine readi(rekord,lancuch,wartosc,default)
480 character*(*) rekord,lancuch
481 integer wartosc,default
484 iread=index(rekord,lancuch)
489 iread=iread+ilen(lancuch)+1
490 read (rekord(iread:),*) wartosc
493 c----------------------------------------------------------------------------
494 subroutine card_concat(card)
496 include 'COMMON.IOUNITS'
498 character*80 karta,ucase
500 read (inp,'(a)') karta
503 do while (karta(80:80).eq.'&')
504 card=card(:ilen(card)+1)//karta(:79)
505 read (inp,'(a)') karta
508 card=card(:ilen(card)+1)//karta
511 c----------------------------------------------------------------------------
520 include 'COMMON.IOUNITS'
521 include 'COMMON.CONTROL'
522 integer lenpre,lenpot,ilen
524 character*16 cformat,cprint
526 integer lenint,lenout
527 call getenv('INPUT',prefix)
528 call getenv('OUTPUT',prefout)
529 call getenv('INTIN',prefintin)
530 call getenv('COORD',cformat)
531 call getenv('PRINTCOOR',cprint)
532 call getenv('SCRATCHDIR',scratchdir)
535 if (index(ucase(cformat),'CX').gt.0) then
542 lenint=ilen(prefintin)
543 C Get the names and open the input files
544 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
546 write (liczba,'(bz,i3.3)') me
547 outname=prefout(:lenout)//'_clust.out_'//liczba
549 outname=prefout(:lenout)//'_clust.out'
552 intinname=prefintin(:lenint)//'.bx'
553 else if (from_cx) then
554 intinname=prefintin(:lenint)//'.cx'
556 intinname=prefintin(:lenint)//'.int'
558 rmsname=prefintin(:lenint)//'.rms'
559 open (jplot,file=prefout(:ilen(prefout))//'.tex',
561 open (jrms,file=rmsname,status='unknown')
562 open(iout,file=outname,status='unknown')
563 C Get parameter filenames and open the parameter files.
564 call getenv('BONDPAR',bondname)
565 open (ibond,file=bondname,status='old')
566 call getenv('THETPAR',thetname)
567 open (ithep,file=thetname,status='old')
568 call getenv('ROTPAR',rotname)
569 open (irotam,file=rotname,status='old')
570 call getenv('TORPAR',torname)
571 open (itorp,file=torname,status='old')
572 call getenv('TORDPAR',tordname)
573 open (itordp,file=tordname,status='old')
574 call getenv('FOURIER',fouriername)
575 open (ifourier,file=fouriername,status='old')
576 call getenv('ELEPAR',elename)
577 open (ielep,file=elename,status='old')
578 call getenv('SIDEPAR',sidename)
579 open (isidep,file=sidename,status='old')
580 call getenv('SIDEP',sidepname)
581 open (isidep1,file=sidepname,status="old")
582 call getenv('SCCORPAR',sccorname)
583 open (isccor,file=sccorname,status="old")
586 C 8/9/01 In the newest version SCp interaction constants are read from a file
587 C Use -DOLDSCP to use hard-coded constants instead.
589 call getenv('SCPPAR',scpname)
590 open (iscpp,file=scpname,status='old')