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
62 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
63 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
64 write (iout,*) "with_dihed_constr ",with_dihed_constr,
65 & " CONSTR_DIST",constr_dist
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'
90 include 'COMMON.TORCNSTR'
94 character*4 sequence(maxres)
95 character*800 weightcard
97 double precision x(maxvar)
98 integer itype_pdb(maxres)
104 C Read weights of the subsequent energy terms.
105 call card_concat(weightcard)
106 call reada(weightcard,'WSC',wsc,1.0d0)
107 call reada(weightcard,'WLONG',wsc,wsc)
108 call reada(weightcard,'WSCP',wscp,1.0d0)
109 call reada(weightcard,'WELEC',welec,1.0D0)
110 call reada(weightcard,'WVDWPP',wvdwpp,welec)
111 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
112 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
113 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
114 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
115 call reada(weightcard,'WTURN3',wturn3,1.0D0)
116 call reada(weightcard,'WTURN4',wturn4,1.0D0)
117 call reada(weightcard,'WTURN6',wturn6,1.0D0)
118 call reada(weightcard,'WSTRAIN',wstrain,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
150 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
151 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
153 10 format (/'Energy-term weights (unscaled):'//
154 & 'WSCC= ',f10.6,' (SC-SC)'/
155 & 'WSCP= ',f10.6,' (SC-p)'/
156 & 'WELEC= ',f10.6,' (p-p electr)'/
157 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
158 & 'WBOND= ',f10.6,' (stretching)'/
159 & 'WANG= ',f10.6,' (bending)'/
160 & 'WSCLOC= ',f10.6,' (SC local)'/
161 & 'WTOR= ',f10.6,' (torsional)'/
162 & 'WTORD= ',f10.6,' (double torsional)'/
163 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
164 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
165 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
166 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
167 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
168 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
169 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
170 & 'WTURN6= ',f10.6,' (turns, 6th order)')
171 if (wcorr4.gt.0.0d0) then
172 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
173 & 'between contact pairs of peptide groups'
174 write (iout,'(2(a,f5.3/))')
175 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
176 & 'Range of quenching the correlation terms:',2*delt_corr
177 else if (wcorr.gt.0.0d0) then
178 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
179 & 'between contact pairs of peptide groups'
181 write (iout,'(a,f8.3)')
182 & 'Scaling factor of 1,4 SC-p interactions:',scal14
183 write (iout,'(a,f8.3)')
184 & 'General scaling factor of SC-p interactions:',scalscp
185 r0_corr=cutoff_corr-delt_corr
187 aad(i,1)=scalscp*aad(i,1)
188 aad(i,2)=scalscp*aad(i,2)
189 bad(i,1)=scalscp*bad(i,1)
190 bad(i,2)=scalscp*bad(i,2)
194 print *,'indpdb=',indpdb,' pdbref=',pdbref
196 C Read sequence if not taken from the pdb file.
197 if (iscode.gt.0) then
198 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
200 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
202 C Convert sequence to numeric code
204 itype(i)=rescode(i,sequence(i),iscode)
207 print '(20i4)',(itype(i),i=1,nres)
211 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
213 if (itype(i).eq.21) then
217 else if (itype(i+1).ne.20) then
219 else if (itype(i).ne.20) then
226 write (iout,*) "ITEL"
228 write (iout,*) i,itype(i),itel(i)
231 print *,'Call Read_Bridge.'
233 if (with_dihed_constr) then
235 read (inp,*) ndih_constr
236 if (ndih_constr.gt.0) then
238 write (iout,*) 'FTORS',ftors
239 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
241 & 'There are',ndih_constr,' constraints on phi angles.'
243 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
246 phi0(i)=deg2rad*phi0(i)
247 drange(i)=deg2rad*drange(i)
254 print *,'NNT=',NNT,' NCT=',NCT
255 if (itype(1).eq.21) nnt=2
256 if (itype(nres).eq.21) nct=nct-1
257 if (nstart.lt.nnt) nstart=nnt
258 if (nend.gt.nct .or. nend.eq.0) nend=nct
259 write (iout,*) "nstart",nstart," nend",nend
262 c read(inp,'(a)') pdbfile
263 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
264 c open(ipdbin,file=pdbfile,status='old',err=33)
266 c 33 write (iout,'(a)') 'Error opening PDB file.'
269 c print *,'Begin reading pdb data'
271 c print *,'Finished reading pdb data'
272 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
274 c itype_pdb(i)=itype(i)
277 c write (iout,'(a,i3)') 'nsup=',nsup
279 c if (nsup.le.(nct-nnt+1)) then
280 c do i=0,nct-nnt+1-nsup
281 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
287 c & 'Error - sequences to be superposed do not match.'
290 c do i=0,nsup-(nct-nnt+1)
291 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
293 c nstart_sup=nstart_sup+i
299 c & 'Error - sequences to be superposed do not match.'
302 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
303 c & ' nstart_seq=',nstart_seq
307 write (iout,*) "molread: REFSTR",refstr
309 if (.not.pdbref) then
310 call read_angles(inp,*38)
312 38 write (iout,'(a)') 'Error reading reference structure.'
314 call mp_stopall(Error_Msg)
316 stop 'Error reading reference structure'
328 call contact(.true.,ncont_ref,icont_ref)
330 c Read distance restraints
331 if (constr_dist.gt.0) then
332 call read_dist_constr
337 c-----------------------------------------------------------------------------
338 logical function seq_comp(itypea,itypeb,length)
340 integer length,itypea(length),itypeb(length)
343 if (itypea(i).ne.itypeb(i)) then
351 c-----------------------------------------------------------------------------
352 subroutine read_bridge
353 C Read information about disulfide bridges.
356 include 'COMMON.IOUNITS'
359 include 'COMMON.INTERACT'
360 include 'COMMON.LOCAL'
361 include 'COMMON.NAMES'
362 include 'COMMON.CHAIN'
363 include 'COMMON.FFIELD'
364 include 'COMMON.SBRIDGE'
365 include 'COMMON.HEADER'
366 include 'COMMON.CONTROL'
367 include 'COMMON.TIME1'
369 include 'COMMON.INFO'
372 C Read bridging residues.
373 read (inp,*) ns,(iss(i),i=1,ns)
375 C Check whether the specified bridging residues are cystines.
377 if (itype(iss(i)).ne.1) then
378 write (iout,'(2a,i3,a)')
379 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
380 & ' can form a disulfide bridge?!!!'
381 write (*,'(2a,i3,a)')
382 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
383 & ' can form a disulfide bridge?!!!'
385 call mp_stopall(error_msg)
391 C Read preformed bridges.
393 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
396 C Check if the residues involved in bridges are in the specified list of
400 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
401 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
402 write (iout,'(a,i3,a)') 'Disulfide pair',i,
403 & ' contains residues present in other pairs.'
404 write (*,'(a,i3,a)') 'Disulfide pair',i,
405 & ' contains residues present in other pairs.'
407 call mp_stopall(error_msg)
414 if (ihpb(i).eq.iss(j)) goto 10
416 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
419 if (jhpb(i).eq.iss(j)) goto 20
421 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
434 c----------------------------------------------------------------------------
435 subroutine read_angles(kanal,*)
440 include 'COMMON.CHAIN'
441 include 'COMMON.IOUNITS'
443 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
444 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
445 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
446 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
448 theta(i)=deg2rad*theta(i)
449 phi(i)=deg2rad*phi(i)
450 alph(i)=deg2rad*alph(i)
451 omeg(i)=deg2rad*omeg(i)
456 c----------------------------------------------------------------------------
457 subroutine reada(rekord,lancuch,wartosc,default)
459 character*(*) rekord,lancuch
460 double precision wartosc,default
463 iread=index(rekord,lancuch)
468 iread=iread+ilen(lancuch)+1
469 read (rekord(iread:),*) wartosc
472 c----------------------------------------------------------------------------
473 subroutine multreada(rekord,lancuch,tablica,dim,default)
476 double precision tablica(dim),default
477 character*(*) rekord,lancuch
483 iread=index(rekord,lancuch)
484 if (iread.eq.0) return
485 iread=iread+ilen(lancuch)+1
486 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
489 c----------------------------------------------------------------------------
490 subroutine readi(rekord,lancuch,wartosc,default)
492 character*(*) rekord,lancuch
493 integer wartosc,default
496 iread=index(rekord,lancuch)
501 iread=iread+ilen(lancuch)+1
502 read (rekord(iread:),*) wartosc
505 c----------------------------------------------------------------------------
506 subroutine multreadi(rekord,lancuch,tablica,dim,default)
509 integer tablica(dim),default
510 character*(*) rekord,lancuch
517 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
518 if (iread.eq.0) return
519 iread=iread+ilen(lancuch)+1
520 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
523 c----------------------------------------------------------------------------
524 subroutine card_concat(card)
526 include 'COMMON.IOUNITS'
528 character*80 karta,ucase
530 read (inp,'(a)') karta
533 do while (karta(80:80).eq.'&')
534 card=card(:ilen(card)+1)//karta(:79)
535 read (inp,'(a)') karta
538 card=card(:ilen(card)+1)//karta
541 c----------------------------------------------------------------------------
550 include 'COMMON.IOUNITS'
551 include 'COMMON.CONTROL'
552 integer lenpre,lenpot,ilen
554 character*16 cformat,cprint
556 integer lenint,lenout
557 call getenv('INPUT',prefix)
558 call getenv('OUTPUT',prefout)
559 call getenv('INTIN',prefintin)
560 call getenv('COORD',cformat)
561 call getenv('PRINTCOOR',cprint)
562 call getenv('SCRATCHDIR',scratchdir)
565 if (index(ucase(cformat),'CX').gt.0) then
572 lenint=ilen(prefintin)
573 C Get the names and open the input files
574 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
576 write (liczba,'(bz,i3.3)') me
577 outname=prefout(:lenout)//'_clust.out_'//liczba
579 outname=prefout(:lenout)//'_clust.out'
582 intinname=prefintin(:lenint)//'.bx'
583 else if (from_cx) then
584 intinname=prefintin(:lenint)//'.cx'
586 intinname=prefintin(:lenint)//'.int'
588 rmsname=prefintin(:lenint)//'.rms'
589 open (jplot,file=prefout(:ilen(prefout))//'.tex',
591 open (jrms,file=rmsname,status='unknown')
592 open(iout,file=outname,status='unknown')
593 C Get parameter filenames and open the parameter files.
594 call getenv('BONDPAR',bondname)
595 open (ibond,file=bondname,status='old')
596 call getenv('THETPAR',thetname)
597 open (ithep,file=thetname,status='old')
598 call getenv('ROTPAR',rotname)
599 open (irotam,file=rotname,status='old')
600 call getenv('TORPAR',torname)
601 open (itorp,file=torname,status='old')
602 call getenv('TORDPAR',tordname)
603 open (itordp,file=tordname,status='old')
604 call getenv('FOURIER',fouriername)
605 open (ifourier,file=fouriername,status='old')
606 call getenv('ELEPAR',elename)
607 open (ielep,file=elename,status='old')
608 call getenv('SIDEPAR',sidename)
609 open (isidep,file=sidename,status='old')
610 call getenv('SIDEP',sidepname)
611 open (isidep1,file=sidepname,status="old")
612 call getenv('SCCORPAR',sccorname)
613 open (isccor,file=sccorname,status="old")
616 C 8/9/01 In the newest version SCp interaction constants are read from a file
617 C Use -DOLDSCP to use hard-coded constants instead.
619 call getenv('SCPPAR',scpname)
620 open (iscpp,file=scpname,status='old')
624 c-------------------------------------------------------------------------------
625 subroutine read_dist_constr
626 implicit real*8 (a-h,o-z)
628 include 'COMMON.CONTROL'
629 include 'COMMON.CHAIN'
630 include 'COMMON.IOUNITS'
631 include 'COMMON.SBRIDGE'
632 integer ifrag_(2,100),ipair_(2,100)
633 double precision wfrag_(100),wpair_(100)
634 character*500 controlcard
635 c write (iout,*) "Calling read_dist_constr"
636 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
637 call card_concat(controlcard)
639 c write (iout,'(a)') controlcard
640 call readi(controlcard,"NFRAG",nfrag_,0)
641 call readi(controlcard,"NPAIR",npair_,0)
642 call readi(controlcard,"NDIST",ndist_,0)
643 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
644 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
645 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
646 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
647 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
648 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
649 write (iout,*) "IFRAG"
651 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
653 write (iout,*) "IPAIR"
655 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
658 if (.not.refstr .and. nfrag_.gt.0) then
660 & "ERROR: no reference structure to compute distance restraints"
662 & "Restraints must be specified explicitly (NDIST=number)"
665 if (nfrag_.lt.2 .and. npair_.gt.0) then
666 write (iout,*) "ERROR: Less than 2 fragments specified",
667 & " but distance restraints between pairs requested"
672 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
673 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
674 & ifrag_(2,i)=nstart_sup+nsup-1
675 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
677 if (wfrag_(i).gt.0.0d0) then
678 do j=ifrag_(1,i),ifrag_(2,i)-1
680 write (iout,*) "j",j," k",k
682 if (constr_dist.eq.1) then
687 forcon(nhpb)=wfrag_(i)
688 else if (constr_dist.eq.2) then
689 if (ddjk.le.dist_cut) then
694 forcon(nhpb)=wfrag_(i)
701 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
703 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
704 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
710 if (wpair_(i).gt.0.0d0) then
718 do j=ifrag_(1,ii),ifrag_(2,ii)
719 do k=ifrag_(1,jj),ifrag_(2,jj)
723 forcon(nhpb)=wpair_(i)
725 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
726 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
732 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
733 & ibecarb(i),forcon(nhpb+1)
734 if (forcon(nhpb+1).gt.0.0d0) then
736 if (ibecarb(i).gt.0) then
740 if (dhpb(nhpb).eq.0.0d0)
741 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
745 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
746 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)