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)
209 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
211 & "Glycine is the first full residue, initial dummy deleted"
217 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
219 & "Glycine is the last full residue, terminal dummy deleted"
223 print '(20i4)',(itype(i),i=1,nres)
227 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
229 if (itype(i).eq.21) then
233 else if (itype(i+1).ne.20) then
235 else if (itype(i).ne.20) then
242 write (iout,*) "ITEL"
244 write (iout,*) i,itype(i),itel(i)
247 print *,'Call Read_Bridge.'
251 print *,'NNT=',NNT,' NCT=',NCT
252 if (itype(1).eq.21) nnt=2
253 if (itype(nres).eq.21) nct=nct-1
254 if (nstart.lt.nnt) nstart=nnt
255 if (nend.gt.nct .or. nend.eq.0) nend=nct
256 write (iout,*) "nstart",nstart," nend",nend
259 c read(inp,'(a)') pdbfile
260 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
261 c open(ipdbin,file=pdbfile,status='old',err=33)
263 c 33 write (iout,'(a)') 'Error opening PDB file.'
266 c print *,'Begin reading pdb data'
268 c print *,'Finished reading pdb data'
269 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
271 c itype_pdb(i)=itype(i)
274 c write (iout,'(a,i3)') 'nsup=',nsup
276 c if (nsup.le.(nct-nnt+1)) then
277 c do i=0,nct-nnt+1-nsup
278 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
284 c & 'Error - sequences to be superposed do not match.'
287 c do i=0,nsup-(nct-nnt+1)
288 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
290 c nstart_sup=nstart_sup+i
296 c & 'Error - sequences to be superposed do not match.'
299 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
300 c & ' nstart_seq=',nstart_seq
304 write (iout,*) "molread: REFSTR",refstr
306 if (.not.pdbref) then
307 call read_angles(inp,*38)
309 38 write (iout,'(a)') 'Error reading reference structure.'
311 call mp_stopall(Error_Msg)
313 stop 'Error reading reference structure'
325 call contact(.true.,ncont_ref,icont_ref)
329 c-----------------------------------------------------------------------------
330 logical function seq_comp(itypea,itypeb,length)
332 integer length,itypea(length),itypeb(length)
335 if (itypea(i).ne.itypeb(i)) then
343 c-----------------------------------------------------------------------------
344 subroutine read_bridge
345 C Read information about disulfide bridges.
348 include 'COMMON.IOUNITS'
351 include 'COMMON.INTERACT'
352 include 'COMMON.LOCAL'
353 include 'COMMON.NAMES'
354 include 'COMMON.CHAIN'
355 include 'COMMON.FFIELD'
356 include 'COMMON.SBRIDGE'
357 include 'COMMON.HEADER'
358 include 'COMMON.CONTROL'
359 include 'COMMON.TIME1'
361 include 'COMMON.INFO'
364 C Read bridging residues.
365 read (inp,*) ns,(iss(i),i=1,ns)
367 C Check whether the specified bridging residues are cystines.
369 if (itype(iss(i)).ne.1) then
370 write (iout,'(2a,i3,a)')
371 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
372 & ' can form a disulfide bridge?!!!'
373 write (*,'(2a,i3,a)')
374 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
375 & ' can form a disulfide bridge?!!!'
377 call mp_stopall(error_msg)
383 C Read preformed bridges.
385 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
388 C Check if the residues involved in bridges are in the specified list of
392 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
393 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
394 write (iout,'(a,i3,a)') 'Disulfide pair',i,
395 & ' contains residues present in other pairs.'
396 write (*,'(a,i3,a)') 'Disulfide pair',i,
397 & ' contains residues present in other pairs.'
399 call mp_stopall(error_msg)
406 if (ihpb(i).eq.iss(j)) goto 10
408 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
411 if (jhpb(i).eq.iss(j)) goto 20
413 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
426 c----------------------------------------------------------------------------
427 subroutine read_angles(kanal,*)
432 include 'COMMON.CHAIN'
433 include 'COMMON.IOUNITS'
435 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
436 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
437 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
438 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
440 theta(i)=deg2rad*theta(i)
441 phi(i)=deg2rad*phi(i)
442 alph(i)=deg2rad*alph(i)
443 omeg(i)=deg2rad*omeg(i)
448 c----------------------------------------------------------------------------
449 subroutine reada(rekord,lancuch,wartosc,default)
451 character*(*) rekord,lancuch
452 double precision wartosc,default
455 iread=index(rekord,lancuch)
460 iread=iread+ilen(lancuch)+1
461 read (rekord(iread:),*) wartosc
464 c----------------------------------------------------------------------------
465 subroutine multreada(rekord,lancuch,tablica,dim,default)
468 double precision tablica(dim),default
469 character*(*) rekord,lancuch
475 iread=index(rekord,lancuch)
476 if (iread.eq.0) return
477 iread=iread+ilen(lancuch)+1
478 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
481 c----------------------------------------------------------------------------
482 subroutine readi(rekord,lancuch,wartosc,default)
484 character*(*) rekord,lancuch
485 integer wartosc,default
488 iread=index(rekord,lancuch)
493 iread=iread+ilen(lancuch)+1
494 read (rekord(iread:),*) wartosc
497 c----------------------------------------------------------------------------
498 subroutine card_concat(card)
500 include 'COMMON.IOUNITS'
502 character*80 karta,ucase
504 read (inp,'(a)') karta
507 do while (karta(80:80).eq.'&')
508 card=card(:ilen(card)+1)//karta(:79)
509 read (inp,'(a)') karta
512 card=card(:ilen(card)+1)//karta
515 c----------------------------------------------------------------------------
524 include 'COMMON.IOUNITS'
525 include 'COMMON.CONTROL'
526 integer lenpre,lenpot,ilen
528 character*16 cformat,cprint
530 integer lenint,lenout
531 call getenv('INPUT',prefix)
532 call getenv('OUTPUT',prefout)
533 call getenv('INTIN',prefintin)
534 call getenv('COORD',cformat)
535 call getenv('PRINTCOOR',cprint)
536 call getenv('SCRATCHDIR',scratchdir)
539 if (index(ucase(cformat),'CX').gt.0) then
546 lenint=ilen(prefintin)
547 C Get the names and open the input files
548 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
550 write (liczba,'(bz,i3.3)') me
551 outname=prefout(:lenout)//'_clust.out_'//liczba
553 outname=prefout(:lenout)//'_clust.out'
556 intinname=prefintin(:lenint)//'.bx'
557 else if (from_cx) then
558 intinname=prefintin(:lenint)//'.cx'
560 intinname=prefintin(:lenint)//'.int'
562 rmsname=prefintin(:lenint)//'.rms'
563 open (jplot,file=prefout(:ilen(prefout))//'.tex',
565 open (jrms,file=rmsname,status='unknown')
566 open(iout,file=outname,status='unknown')
567 C Get parameter filenames and open the parameter files.
568 call getenv('BONDPAR',bondname)
569 open (ibond,file=bondname,status='old')
570 call getenv('THETPAR',thetname)
571 open (ithep,file=thetname,status='old')
572 call getenv('ROTPAR',rotname)
573 open (irotam,file=rotname,status='old')
574 call getenv('TORPAR',torname)
575 open (itorp,file=torname,status='old')
576 call getenv('TORDPAR',tordname)
577 open (itordp,file=tordname,status='old')
578 call getenv('FOURIER',fouriername)
579 open (ifourier,file=fouriername,status='old')
580 call getenv('ELEPAR',elename)
581 open (ielep,file=elename,status='old')
582 call getenv('SIDEPAR',sidename)
583 open (isidep,file=sidename,status='old')
584 call getenv('SIDEP',sidepname)
585 open (isidep1,file=sidepname,status="old")
586 call getenv('SCCORPAR',sccorname)
587 open (isccor,file=sccorname,status="old")
590 C 8/9/01 In the newest version SCp interaction constants are read from a file
591 C Use -DOLDSCP to use hard-coded constants instead.
593 call getenv('SCPPAR',scpname)
594 open (iscpp,file=scpname,status='old')