subroutine read_control C C Read molecular data C implicit none include 'DIMENSIONS' include 'sizesclu.dat' include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.HEADER' include 'COMMON.FFIELD' include 'COMMON.FREE' include 'COMMON.INTERACT' include "COMMON.SPLITELE" character*320 controlcard,ucase #ifdef MPL include 'COMMON.INFO' #endif integer i read (INP,'(a80)') titel call card_concat(controlcard) call readi(controlcard,'NRES',nres,0) call readi(controlcard,'RESCALE',rescale_mode,2) call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0) write (iout,*) "DISTCHAINMAX",distchainmax C Reading the dimensions of box in x,y,z coordinates call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) call readi(controlcard,'PDBOUT',outpdb,0) call readi(controlcard,'MOL2OUT',outmol2,0) refstr=(index(controlcard,'REFSTR').gt.0) write (iout,*) "REFSTR",refstr pdbref=(index(controlcard,'PDBREF').gt.0) iscode=index(controlcard,'ONE_LETTER') tree=(index(controlcard,'MAKE_TREE').gt.0) min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) call readi(controlcard,'NCUT',ncut,1) call readi(controlcard,'SYM',symetr,1) write (iout,*) 'sym', symetr call readi(controlcard,'NSTART',nstart,0) call readi(controlcard,'NEND',nend,0) call reada(controlcard,'ECUT',ecut,10.0d0) call reada(controlcard,'PROB',prob_limit,0.99d0) write (iout,*) "Probability limit",prob_limit lgrp=(index(controlcard,'LGRP').gt.0) caonly=(index(controlcard,'CA_ONLY').gt.0) print_dist=(index(controlcard,'PRINT_DIST').gt.0) call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) call readi(controlcard,'IOPT',iopt,2) lside = index(controlcard,"SIDE").gt.0 efree = index(controlcard,"EFREE").gt.0 call readi(controlcard,'NTEMP',nT,1) write (iout,*) "nT",nT call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0) write (iout,*) "nT",nT write (iout,*) 'beta_h',(beta_h(i),i=1,nT) do i=1,nT beta_h(i)=1.0d0/(1.987D-3*beta_h(i)) enddo write (iout,*) 'beta_h',(beta_h(i),i=1,nT) lprint_cart=index(controlcard,"PRINT_CART") .gt.0 lprint_int=index(controlcard,"PRINT_INT") .gt.0 if (min_var) iopt=1 return end c-------------------------------------------------------------------------- subroutine molread C C Read molecular data. C implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.CONTACTS' include 'COMMON.TIME1' #ifdef MPL include 'COMMON.INFO' #endif character*4 sequence(maxres) character*800 weightcard integer rescode double precision x(maxvar) integer itype_pdb(maxres) logical seq_comp integer i,j,kkk C C Body C C Read weights of the subsequent energy terms. call card_concat(weightcard) call reada(weightcard,'WSC',wsc,1.0d0) call reada(weightcard,'WLONG',wsc,wsc) call reada(weightcard,'WSCP',wscp,1.0d0) call reada(weightcard,'WELEC',welec,1.0D0) call reada(weightcard,'WVDWPP',wvdwpp,welec) call reada(weightcard,'WEL_LOC',wel_loc,1.0D0) call reada(weightcard,'WCORR4',wcorr4,0.0D0) call reada(weightcard,'WCORR5',wcorr5,0.0D0) call reada(weightcard,'WCORR6',wcorr6,0.0D0) call reada(weightcard,'WTURN3',wturn3,1.0D0) call reada(weightcard,'WTURN4',wturn4,1.0D0) call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 weights(1)=wsc weights(2)=wscp weights(3)=welec weights(4)=wcorr weights(5)=wcorr5 weights(6)=wcorr6 weights(7)=wel_loc weights(8)=wturn3 weights(9)=wturn4 weights(10)=wturn6 weights(11)=wang weights(12)=wscloc weights(13)=wtor weights(14)=wtor_d weights(15)=wstrain weights(16)=wvdwpp weights(17)=wbond weights(18)=scal14 weights(19)=wsccor write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3, & wturn4,wturn6,wsccor 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ & 'WELEC= ',f10.6,' (p-p electr)'/ & 'WVDWPP= ',f10.6,' (p-p VDW)'/ & 'WBOND= ',f10.6,' (stretching)'/ & 'WANG= ',f10.6,' (bending)'/ & 'WSCLOC= ',f10.6,' (SC local)'/ & 'WTOR= ',f10.6,' (torsional)'/ & 'WTORD= ',f10.6,' (double torsional)'/ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ & 'WCORR4= ',f10.6,' (multi-body 4th order)'/ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ & 'WTURN6= ',f10.6,' (turns, 6th order)'/ & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)') if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', & 'between contact pairs of peptide groups' write (iout,'(2(a,f5.3/))') & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr, & 'Range of quenching the correlation terms:',2*delt_corr else if (wcorr.gt.0.0d0) then write (iout,'(/2a/)') 'Hydrogen-bonding correlation ', & 'between contact pairs of peptide groups' endif write (iout,'(a,f8.3)') & 'Scaling factor of 1,4 SC-p interactions:',scal14 write (iout,'(a,f8.3)') & 'General scaling factor of SC-p interactions:',scalscp r0_corr=cutoff_corr-delt_corr do i=1,20 aad(i,1)=scalscp*aad(i,1) aad(i,2)=scalscp*aad(i,2) bad(i,1)=scalscp*bad(i,1) bad(i,2)=scalscp*bad(i,2) enddo call flush(iout) print *,'indpdb=',indpdb,' pdbref=',pdbref C Read sequence if not taken from the pdb file. if (iscode.gt.0) then read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres) else read (inp,'(20(1x,a3))') (sequence(i),i=1,nres) endif C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo print *,nres print '(20i4)',(itype(i),i=1,nres) do i=1,nres #ifdef PROCOR if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then #else if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR else if (iabs(itype(i+1)).ne.20) then #else else if (iabs(itype(i)).ne.20) then #endif itel(i)=1 else itel(i)=2 endif enddo write (iout,*) "ITEL" do i=1,nres-1 write (iout,*) i,itype(i),itel(i) enddo print *,'Call Read_Bridge.' call read_bridge nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT if (itype(1).eq.ntyp1) nnt=2 if (itype(nres).eq.ntyp1) nct=nct-1 if (nstart.lt.nnt) nstart=nnt if (nend.gt.nct .or. nend.eq.0) nend=nct write (iout,*) "nstart",nstart," nend",nend nres0=nres c if (pdbref) then c read(inp,'(a)') pdbfile c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile c open(ipdbin,file=pdbfile,status='old',err=33) c goto 34 c 33 write (iout,'(a)') 'Error opening PDB file.' c stop c 34 continue c print *,'Begin reading pdb data' c call readpdb c print *,'Finished reading pdb data' c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup c do i=1,nres c itype_pdb(i)=itype(i) c enddo c close (ipdbin) c write (iout,'(a,i3)') 'nsup=',nsup c nstart_seq=nnt c if (nsup.le.(nct-nnt+1)) then c do i=0,nct-nnt+1-nsup c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then c nstart_seq=nnt+i c goto 111 c endif c enddo c write (iout,'(a)') c & 'Error - sequences to be superposed do not match.' c stop c else c do i=0,nsup-(nct-nnt+1) c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) c & then c nstart_sup=nstart_sup+i c nsup=nct-nnt+1 c goto 111 c endif c enddo c write (iout,'(a)') c & 'Error - sequences to be superposed do not match.' c endif c 111 continue c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup, c & ' nstart_seq=',nstart_seq c endif call init_int_table call setup_var write (iout,*) "molread: REFSTR",refstr if (refstr) then if (.not.pdbref) then call read_angles(inp,*38) goto 39 38 write (iout,'(a)') 'Error reading reference structure.' #ifdef MPL call mp_stopall(Error_Msg) #else stop 'Error reading reference structure' #endif 39 call chainbuild nstart_sup=nnt nstart_seq=nnt nsup=nct-nnt+1 kkk=1 do i=1,2*nres do j=1,3 cref(j,i,kkk)=c(j,i) enddo enddo endif call contact(.true.,ncont_ref,icont_ref) endif return end c----------------------------------------------------------------------------- logical function seq_comp(itypea,itypeb,length) implicit none integer length,itypea(length),itypeb(length) integer i do i=1,length if (itypea(i).ne.itypeb(i)) then seq_comp=.false. return endif enddo seq_comp=.true. return end c----------------------------------------------------------------------------- subroutine read_bridge C Read information about disulfide bridges. implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.TIME1' #ifdef MPL include 'COMMON.INFO' #endif integer i,j C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) print *,'ns=',ns C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then write (iout,'(2a,i3,a)') & 'Do you REALLY think that the residue ',restyp(iss(i)),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') & 'Do you REALLY think that the residue ',restyp(iss(i)),i, & ' can form a disulfide bridge?!!!' #ifdef MPL call mp_stopall(error_msg) #else stop #endif endif enddo C Read preformed bridges. if (ns.gt.0) then read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) if (nss.gt.0) then nhpb=nss C Check if the residues involved in bridges are in the specified list of C bridging residues. do i=1,nss do j=1,i-1 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then write (iout,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' write (*,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' #ifdef MPL call mp_stopall(error_msg) #else stop #endif endif enddo do j=1,ns if (ihpb(i).eq.iss(j)) goto 10 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 10 continue do j=1,ns if (jhpb(i).eq.iss(j)) goto 20 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 20 continue dhpb(i)=dbr forcon(i)=fbr enddo do i=1,nss ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres enddo endif endif return end c---------------------------------------------------------------------------- subroutine read_angles(kanal,*) implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' integer i,kanal read (kanal,*,err=10,end=10) (theta(i),i=3,nres) read (kanal,*,err=10,end=10) (phi(i),i=4,nres) read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1) read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1) do i=1,nres theta(i)=deg2rad*theta(i) phi(i)=deg2rad*phi(i) alph(i)=deg2rad*alph(i) omeg(i)=deg2rad*omeg(i) enddo return 10 return1 end c---------------------------------------------------------------------------- subroutine reada(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch double precision wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*) wartosc return end c---------------------------------------------------------------------------- subroutine multreada(rekord,lancuch,tablica,dim,default) implicit none integer dim,i double precision tablica(dim),default character*(*) rekord,lancuch integer ilen,iread external ilen do i=1,dim tablica(i)=default enddo iread=index(rekord,lancuch) if (iread.eq.0) return iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end c---------------------------------------------------------------------------- subroutine readi(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch integer wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*) wartosc return end c---------------------------------------------------------------------------- subroutine card_concat(card) include 'DIMENSIONS' include 'COMMON.IOUNITS' character*(*) card character*80 karta,ucase external ilen read (inp,'(a)') karta karta=ucase(karta) card=' ' do while (karta(80:80).eq.'&') card=card(:ilen(card)+1)//karta(:79) read (inp,'(a)') karta karta=ucase(karta) enddo card=card(:ilen(card)+1)//karta return end c---------------------------------------------------------------------------- subroutine openunits implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" character*3 liczba include "COMMON.MPI" #endif include 'COMMON.IOUNITS' include 'COMMON.CONTROL' integer lenpre,lenpot,ilen external ilen character*16 cformat,cprint character*16 ucase integer lenint,lenout call getenv('INPUT',prefix) call getenv('OUTPUT',prefout) call getenv('INTIN',prefintin) call getenv('COORD',cformat) call getenv('PRINTCOOR',cprint) call getenv('SCRATCHDIR',scratchdir) from_bx=.true. from_cx=.false. if (index(ucase(cformat),'CX').gt.0) then from_cx=.true. from_bx=.false. endif from_cart=.true. lenpre=ilen(prefix) lenout=ilen(prefout) lenint=ilen(prefintin) C Get the names and open the input files open (inp,file=prefix(:ilen(prefix))//'.inp',status='old') #ifdef MPI write (liczba,'(bz,i3.3)') me outname=prefout(:lenout)//'_clust.out_'//liczba #else outname=prefout(:lenout)//'_clust.out' #endif if (from_bx) then intinname=prefintin(:lenint)//'.bx' else if (from_cx) then intinname=prefintin(:lenint)//'.cx' else intinname=prefintin(:lenint)//'.int' endif rmsname=prefintin(:lenint)//'.rms' open (jplot,file=prefout(:ilen(prefout))//'.tex', & status='unknown') open (jrms,file=rmsname,status='unknown') open(iout,file=outname,status='unknown') C Get parameter filenames and open the parameter files. call getenv('BONDPAR',bondname) open (ibond,file=bondname,status='old') call getenv('THETPAR',thetname) open (ithep,file=thetname,status='old') call getenv('ROTPAR',rotname) open (irotam,file=rotname,status='old') call getenv('TORPAR',torname) open (itorp,file=torname,status='old') call getenv('TORDPAR',tordname) open (itordp,file=tordname,status='old') call getenv('FOURIER',fouriername) open (ifourier,file=fouriername,status='old') call getenv('ELEPAR',elename) open (ielep,file=elename,status='old') call getenv('SIDEPAR',sidename) open (isidep,file=sidename,status='old') call getenv('SIDEP',sidepname) open (isidep1,file=sidepname,status="old") call getenv('SCCORPAR',sccorname) open (isccor,file=sccorname,status="old") #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file C Use -DOLDSCP to use hard-coded constants instead. C call getenv('SCPPAR',scpname) open (iscpp,file=scpname,status='old') #endif return end