subroutine read_control C C Read molecular data C include 'DIMENSIONS' include 'sizesclu.dat' include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.THREAD' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.HEADER' character*320 controlcard,ucase #ifdef MPL include 'COMMON.INFO' #endif read (INP,'(a80)') titel call card_concat(controlcard) C Set up the time limit (caution! The time must be input in minutes!) call reada(controlcard,'TIMLIM',timlim,60.0D0) write (iout,'(a,f10.1)') 'Time limit (min):',timlim timlim=60.0D0*timlim #ifdef MPL timlim=timlim/nprocs #endif modecalc=0 call readi(controlcard,'NRES',nres,0) call reada(controlcard,'TEMP',betaT,298.0d0) write (iout,*) "temperature",betaT betaT=1.0d0/(betaT*1.987D-3) write (iout,*) "betaT",betaT call readi(controlcard,'PDBOUT',outpdb,0) call readi(controlcard,'MOL2OUT',outmol2,0) refstr=(index(controlcard,'REFSTR').gt.0) 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 reada(controlcard,'ECUT',ecut,10.0d0) 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 if (min_var) iopt=1 return end c-------------------------------------------------------------------------- subroutine molread C C Read molecular data. C 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.DBASE' include 'COMMON.THREAD' include 'COMMON.CONTACTS' include 'COMMON.TIME1' #ifdef MPL include 'COMMON.INFO' #endif character*4 sequence(maxres) integer rescode double precision x(maxvar) character*64 pdbfile dimension itype_pdb(maxres) logical seq_comp C C Body C if (indpdb.gt.0 .or. pdbref) then read(inp,'(a)') pdbfile write (iout,'(2a)') 'PDB data will be read from file ',pdbfile open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a)') 'Error opening PDB file.' stop 34 continue print *,'Begin reading pdb data' call readpdb print *,'Finished reading pdb data' write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup do i=1,nres itype_pdb(i)=itype(i) enddo close (ipdbin) endif if (indpdb.eq.0) then 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 endif print *,nres print '(20i4)',(itype(i),i=1,nres) do i=1,nres if (itype(i).ne.20) then itel(i)=1 else itel(i)=2 endif enddo print *,'Call Read_Bridge.' call read_bridge nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT if (itype(1).eq.21) nnt=2 if (itype(nres).eq.21) nct=nct-1 if (pdbref) then write (iout,'(a,i3)') 'nsup=',nsup nstart_seq=nnt if (nsup.le.(nct-nnt+1)) then do i=0,nct-nnt+1-nsup if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then nstart_seq=nnt+i goto 111 endif enddo write (iout,'(a)') & 'Error - sequences to be superposed do not match.' stop else do i=0,nsup-(nct-nnt+1) if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) & then nstart_sup=nstart_sup+i nsup=nct-nnt+1 goto 111 endif enddo write (iout,'(a)') & 'Error - sequences to be superposed do not match.' endif 111 continue write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup, & ' nstart_seq=',nstart_seq endif 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 do i=1,2*nres do j=1,3 cref(j,i)=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. 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.DBASE' include 'COMMON.THREAD' include 'COMMON.TIME1' #ifdef MPL include 'COMMON.INFO' #endif 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,*) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' 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 include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' integer lenpre,lenpot,ilen external ilen character*16 cfrom_pdb,cprint character*16 ucase call getenv('INPUT',prefix) call getenv('OUTPUT',prefout) call getenv('INTIN',prefintin) call getenv('PDB',cfrom_pdb) call getenv('PRINTCOOR',cprint) from_cart = index(ucase(cfrom_pdb),'CART').gt.0 lprint_cart = index(ucase(cprint),'PRINT_CART').gt.0 lprint_int = index(ucase(cprint),'PRINT_INT').gt.0 if (from_cart .and. .not.lprint_int) then lprint_cart=.true. lprint_int=.false. endif if (.not.lprint_cart .and. .not.lprint_int) lprint_int=.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') outname=prefout(:lenout)//'_clust.out' if (lprint_cart) then intname=prefintin(:lenint)//'_clust'//'.x' intinname=prefintin(:lenint)//'.x' else intname=prefintin(:lenint)//'_clust'//'.int' intinname=prefintin(:lenint)//'.int' endif rmsname=prefintin(:lenint)//'.rms' statinname=prefintin(:lenint)//'.stat' print *,statinname statoutname=prefintin(:lenint)//'_clust.stat' open (jplot,file=prefout(:ilen(prefout))//'.tex', & status='unknown') print *,'unit',jplot,' opened' open (intin,file=intinname,status='old') print *,'unit',intin,' opened' open (jrms,file=rmsname,status='unknown') open (jstatin,file=statinname,status='unknown') open (jstatout,file=statoutname,status='unknown') #ifdef AIX open(iout,file=outname,status='unknown') print *,'unit',iout,' opened' open(igeom,file=intname,status='unknown',position='append') print *,'unit',igeom,' opened' #elif defined(G77) open(iout,file=outname,status='unknown') open(igeom,file=intname,status='unknown',access='append') #else open(iout,file=outname,status='unknown') open(igeom,file=intname,status='unknown',position='append') #endif return end