source/cluster/unres/src/readrtns.F

   1       subroutine read_control
   2 C
   3 C Read molecular data
   4 C
   5       include 'DIMENSIONS'
   6       include 'sizesclu.dat'
   7       include 'COMMON.IOUNITS'
   8       include 'COMMON.TIME1'
   9       include 'COMMON.THREAD'
  10       include 'COMMON.SBRIDGE'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.HEADER'
  15       character*320 controlcard,ucase
  16 #ifdef MPL
  17       include 'COMMON.INFO'
  18 #endif
  19
  20       read (INP,'(a80)') titel
  21       call card_concat(controlcard)
  22
  23 C Set up the time limit (caution! The time must be input in minutes!)
  24       call reada(controlcard,'TIMLIM',timlim,60.0D0)
  25       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
  26       timlim=60.0D0*timlim
  27 #ifdef MPL
  28       timlim=timlim/nprocs
  29 #endif
  30
  31       modecalc=0
  32       call readi(controlcard,'NRES',nres,0)
  33       call reada(controlcard,'TEMP',betaT,298.0d0)
  34       write (iout,*) "temperature",betaT
  35       betaT=1.0d0/(betaT*1.987D-3)
  36       write (iout,*) "betaT",betaT
  37       call readi(controlcard,'PDBOUT',outpdb,0)
  38       call readi(controlcard,'MOL2OUT',outmol2,0)
  39       refstr=(index(controlcard,'REFSTR').gt.0)
  40       pdbref=(index(controlcard,'PDBREF').gt.0)
  41       iscode=index(controlcard,'ONE_LETTER')
  42       tree=(index(controlcard,'MAKE_TREE').gt.0)
  43       min_var=(index(controlcard,'MINVAR').gt.0)
  44       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
  45       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
  46       call readi(controlcard,'NCUT',ncut,1)
  47       call reada(controlcard,'ECUT',ecut,10.0d0)
  48       lgrp=(index(controlcard,'LGRP').gt.0)
  49       caonly=(index(controlcard,'CA_ONLY').gt.0)
  50       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
  51       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
  52       call readi(controlcard,'IOPT',iopt,2)
  53       lside = index(controlcard,"SIDE").gt.0
  54       efree = index(controlcard,"EFREE").gt.0
  55       if (min_var) iopt=1
  56       return
  57       end
  58 c--------------------------------------------------------------------------
  59       subroutine molread
  60 C
  61 C Read molecular data.
  62 C
  63       include 'DIMENSIONS'
  64       include 'COMMON.IOUNITS'
  65       include 'COMMON.GEO'
  66       include 'COMMON.VAR'
  67       include 'COMMON.INTERACT'
  68       include 'COMMON.LOCAL'
  69       include 'COMMON.NAMES'
  70       include 'COMMON.CHAIN'
  71       include 'COMMON.FFIELD'
  72       include 'COMMON.SBRIDGE'
  73       include 'COMMON.HEADER'
  74       include 'COMMON.CONTROL'
  75       include 'COMMON.DBASE'
  76       include 'COMMON.THREAD'
  77       include 'COMMON.CONTACTS'
  78       include 'COMMON.TIME1'
  79 #ifdef MPL
  80       include 'COMMON.INFO'
  81 #endif
  82       character*4 sequence(maxres)
  83       integer rescode
  84       double precision x(maxvar)
  85       character*64 pdbfile
  86       dimension itype_pdb(maxres)
  87       logical seq_comp
  88 C
  89 C Body
  90 C
  91       if (indpdb.gt.0 .or. pdbref) then
  92         read(inp,'(a)') pdbfile
  93         write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
  94         open(ipdbin,file=pdbfile,status='old',err=33)
  95         goto 34
  96   33    write (iout,'(a)') 'Error opening PDB file.'
  97         stop
  98   34    continue
  99         print *,'Begin reading pdb data'
 100         call readpdb
 101         print *,'Finished reading pdb data'
 102         write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 103         do i=1,nres
 104           itype_pdb(i)=itype(i)
 105         enddo
 106         close (ipdbin)
 107       endif
 108       if (indpdb.eq.0) then
 109 C Read sequence if not taken from the pdb file.
 110         if (iscode.gt.0) then
 111           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 112         else
 113           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 114         endif
 115 C Convert sequence to numeric code
 116         do i=1,nres
 117           itype(i)=rescode(i,sequence(i),iscode)
 118         enddo
 119       endif
 120       print *,nres
 121       print '(20i4)',(itype(i),i=1,nres)
 122       do i=1,nres
 123         if (itype(i).ne.20) then
 124           itel(i)=1
 125         else
 126           itel(i)=2
 127         endif
 128       enddo
 129       print *,'Call Read_Bridge.'
 130       call read_bridge
 131       nnt=1
 132       nct=nres
 133       print *,'NNT=',NNT,' NCT=',NCT
 134       if (itype(1).eq.21) nnt=2
 135       if (itype(nres).eq.21) nct=nct-1
 136       if (pdbref) then
 137         write (iout,'(a,i3)') 'nsup=',nsup
 138         nstart_seq=nnt
 139         if (nsup.le.(nct-nnt+1)) then
 140           do i=0,nct-nnt+1-nsup
 141             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 142               nstart_seq=nnt+i
 143               goto 111
 144             endif
 145           enddo
 146           write (iout,'(a)')
 147      &            'Error - sequences to be superposed do not match.'
 148           stop
 149         else
 150           do i=0,nsup-(nct-nnt+1)
 151             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 152      &      then
 153               nstart_sup=nstart_sup+i
 154               nsup=nct-nnt+1
 155               goto 111
 156             endif
 157           enddo
 158           write (iout,'(a)')
 159      &            'Error - sequences to be superposed do not match.'
 160         endif
 161   111   continue
 162         write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 163      &                 ' nstart_seq=',nstart_seq
 164       endif
 165       if (refstr) then
 166         if (.not.pdbref) then
 167           call read_angles(inp,*38)
 168           goto 39
 169    38     write (iout,'(a)') 'Error reading reference structure.'
 170 #ifdef MPL
 171           call mp_stopall(Error_Msg)
 172 #else
 173           stop 'Error reading reference structure'
 174 #endif
 175    39     call chainbuild
 176           nstart_sup=nnt
 177           nstart_seq=nnt
 178           nsup=nct-nnt+1
 179           do i=1,2*nres
 180             do j=1,3
 181               cref(j,i)=c(j,i)
 182             enddo
 183           enddo
 184         endif
 185         call contact(.true.,ncont_ref,icont_ref)
 186       endif
 187       return
 188       end
 189 c-----------------------------------------------------------------------------
 190       logical function seq_comp(itypea,itypeb,length)
 191       implicit none
 192       integer length,itypea(length),itypeb(length)
 193       integer i
 194       do i=1,length
 195         if (itypea(i).ne.itypeb(i)) then
 196           seq_comp=.false.
 197           return
 198         endif
 199       enddo
 200       seq_comp=.true.
 201       return
 202       end
 203 c-----------------------------------------------------------------------------
 204       subroutine read_bridge
 205 C Read information about disulfide bridges.
 206       include 'DIMENSIONS'
 207       include 'COMMON.IOUNITS'
 208       include 'COMMON.GEO'
 209       include 'COMMON.VAR'
 210       include 'COMMON.INTERACT'
 211       include 'COMMON.LOCAL'
 212       include 'COMMON.NAMES'
 213       include 'COMMON.CHAIN'
 214       include 'COMMON.FFIELD'
 215       include 'COMMON.SBRIDGE'
 216       include 'COMMON.HEADER'
 217       include 'COMMON.CONTROL'
 218       include 'COMMON.DBASE'
 219       include 'COMMON.THREAD'
 220       include 'COMMON.TIME1'
 221 #ifdef MPL
 222       include 'COMMON.INFO'
 223 #endif
 224 C Read bridging residues.
 225       read (inp,*) ns,(iss(i),i=1,ns)
 226       print *,'ns=',ns
 227 C Check whether the specified bridging residues are cystines.
 228       do i=1,ns
 229         if (itype(iss(i)).ne.1) then
 230           write (iout,'(2a,i3,a)')
 231      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 232      &   ' can form a disulfide bridge?!!!'
 233           write (*,'(2a,i3,a)')
 234      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 235      &   ' can form a disulfide bridge?!!!'
 236 #ifdef MPL
 237          call mp_stopall(error_msg)
 238 #else
 239          stop
 240 #endif
 241         endif
 242       enddo
 243 C Read preformed bridges.
 244       if (ns.gt.0) then
 245       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 246       if (nss.gt.0) then
 247         nhpb=nss
 248 C Check if the residues involved in bridges are in the specified list of
 249 C bridging residues.
 250         do i=1,nss
 251           do j=1,i-1
 252             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 253      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 254               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 255      &      ' contains residues present in other pairs.'
 256               write (*,'(a,i3,a)') 'Disulfide pair',i,
 257      &      ' contains residues present in other pairs.'
 258 #ifdef MPL
 259               call mp_stopall(error_msg)
 260 #else
 261               stop
 262 #endif
 263             endif
 264           enddo
 265           do j=1,ns
 266             if (ihpb(i).eq.iss(j)) goto 10
 267           enddo
 268           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 269    10     continue
 270           do j=1,ns
 271             if (jhpb(i).eq.iss(j)) goto 20
 272           enddo
 273           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 274    20     continue
 275           dhpb(i)=dbr
 276           forcon(i)=fbr
 277         enddo
 278         do i=1,nss
 279           ihpb(i)=ihpb(i)+nres
 280           jhpb(i)=jhpb(i)+nres
 281         enddo
 282       endif
 283       endif
 284       return
 285       end
 286 c----------------------------------------------------------------------------
 287       subroutine read_angles(kanal,*)
 288       include 'DIMENSIONS'
 289       include 'COMMON.GEO'
 290       include 'COMMON.VAR'
 291       include 'COMMON.CHAIN'
 292       include 'COMMON.IOUNITS'
 293       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 294       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 295       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 296       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 297       do i=1,nres
 298         theta(i)=deg2rad*theta(i)
 299         phi(i)=deg2rad*phi(i)
 300         alph(i)=deg2rad*alph(i)
 301         omeg(i)=deg2rad*omeg(i)
 302       enddo
 303       return
 304    10 return1
 305       end
 306 c----------------------------------------------------------------------------
 307       subroutine reada(rekord,lancuch,wartosc,default)
 308       implicit none
 309       character*(*) rekord,lancuch
 310       double precision wartosc,default
 311       integer ilen,iread
 312       external ilen
 313       iread=index(rekord,lancuch)
 314       if (iread.eq.0) then
 315         wartosc=default
 316         return
 317       endif
 318       iread=iread+ilen(lancuch)+1
 319       read (rekord(iread:),*) wartosc
 320       return
 321       end
 322 c----------------------------------------------------------------------------
 323       subroutine multreada(rekord,lancuch,tablica,dim,default)
 324       implicit none
 325       integer dim,i
 326       double precision tablica(dim),default
 327       character*(*) rekord,lancuch
 328       integer ilen,iread
 329       external ilen
 330       do i=1,dim
 331         tablica(i)=default
 332       enddo
 333       iread=index(rekord,lancuch)
 334       if (iread.eq.0) return
 335       iread=iread+ilen(lancuch)+1
 336       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 337    10 return
 338       end
 339 c----------------------------------------------------------------------------
 340       subroutine readi(rekord,lancuch,wartosc,default)
 341       implicit none
 342       character*(*) rekord,lancuch
 343       integer wartosc,default
 344       integer ilen,iread
 345       external ilen
 346       iread=index(rekord,lancuch)
 347       if (iread.eq.0) then
 348         wartosc=default
 349         return
 350       endif
 351       iread=iread+ilen(lancuch)+1
 352       read (rekord(iread:),*) wartosc
 353       return
 354       end
 355 c----------------------------------------------------------------------------
 356       subroutine card_concat(card)
 357       include 'DIMENSIONS'
 358       include 'COMMON.IOUNITS'
 359       character*(*) card
 360       character*80 karta,ucase
 361       external ilen
 362       read (inp,'(a)') karta
 363       karta=ucase(karta)
 364       card=' '
 365       do while (karta(80:80).eq.'&')
 366         card=card(:ilen(card)+1)//karta(:79)
 367         read (inp,'(a)') karta
 368         karta=ucase(karta)
 369       enddo
 370       card=card(:ilen(card)+1)//karta
 371       return
 372       end
 373 c----------------------------------------------------------------------------
 374       subroutine openunits
 375       include 'DIMENSIONS'
 376       include 'COMMON.IOUNITS'
 377       include 'COMMON.CONTROL'
 378       integer lenpre,lenpot,ilen
 379       external ilen
 380       character*16 cfrom_pdb,cprint
 381       character*16 ucase
 382       call getenv('INPUT',prefix)
 383       call getenv('OUTPUT',prefout)
 384       call getenv('INTIN',prefintin)
 385       call getenv('PDB',cfrom_pdb)
 386       call getenv('PRINTCOOR',cprint)
 387       from_cart = index(ucase(cfrom_pdb),'CART').gt.0
 388       lprint_cart = index(ucase(cprint),'PRINT_CART').gt.0
 389       lprint_int = index(ucase(cprint),'PRINT_INT').gt.0
 390       if (from_cart  .and. .not.lprint_int) then
 391         lprint_cart=.true.
 392         lprint_int=.false.
 393       endif
 394       if (.not.lprint_cart .and. .not.lprint_int) lprint_int=.true.
 395       lenpre=ilen(prefix)
 396       lenout=ilen(prefout)
 397       lenint=ilen(prefintin)
 398 C Get the names and open the input files
 399       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 400       outname=prefout(:lenout)//'_clust.out'
 401       if (lprint_cart) then
 402         intname=prefintin(:lenint)//'_clust'//'.x'
 403         intinname=prefintin(:lenint)//'.x'
 404       else
 405         intname=prefintin(:lenint)//'_clust'//'.int'
 406         intinname=prefintin(:lenint)//'.int'
 407       endif
 408       rmsname=prefintin(:lenint)//'.rms'
 409       statinname=prefintin(:lenint)//'.stat'
 410       print *,statinname
 411       statoutname=prefintin(:lenint)//'_clust.stat'
 412       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 413      &       status='unknown')
 414       print *,'unit',jplot,' opened'
 415       open (intin,file=intinname,status='old')
 416       print *,'unit',intin,' opened'
 417       open (jrms,file=rmsname,status='unknown')
 418       open (jstatin,file=statinname,status='unknown')
 419       open (jstatout,file=statoutname,status='unknown')
 420 #ifdef AIX
 421       open(iout,file=outname,status='unknown')
 422       print *,'unit',iout,' opened'
 423       open(igeom,file=intname,status='unknown',position='append')
 424       print *,'unit',igeom,' opened'
 425 #elif defined(G77)
 426       open(iout,file=outname,status='unknown')
 427       open(igeom,file=intname,status='unknown',access='append')
 428 #else
 429       open(iout,file=outname,status='unknown')
 430       open(igeom,file=intname,status='unknown',position='append')
 431 #endif
 432       return
 433       end