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       print *,"print_dist",print_dist
  52       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
  53       call readi(controlcard,'IOPT',iopt,2)
  54       print *,"iopt",iopt
  55       lside = index(controlcard,"SIDE").gt.0
  56       efree = index(controlcard,"EFREE").gt.0
  57       print *,controlcard
  58       call readi(controlcard,'ISTART',is,1)
  59       print *,"is",is
  60       call readi(controlcard,'IEND',ie,10000000)
  61       print *,"ie",ie
  62       call readi(controlcard,'ISAMPL',isampl,1)
  63       print *,"isampl",isampl
  64       call reada(controlcard,'TS',ts,1.0d10)
  65       call reada(controlcard,'TE',te,1.0d10)
  66       print *,"is",is," ie",ie
  67       if (min_var) iopt=1
  68       return
  69       end
  70 c--------------------------------------------------------------------------
  71       subroutine molread
  72 C
  73 C Read molecular data.
  74 C
  75       include 'DIMENSIONS'
  76       include 'COMMON.IOUNITS'
  77       include 'COMMON.GEO'
  78       include 'COMMON.VAR'
  79       include 'COMMON.INTERACT'
  80       include 'COMMON.LOCAL'
  81       include 'COMMON.NAMES'
  82       include 'COMMON.CHAIN'
  83       include 'COMMON.FFIELD'
  84       include 'COMMON.SBRIDGE'
  85       include 'COMMON.HEADER'
  86       include 'COMMON.CONTROL'
  87       include 'COMMON.DBASE'
  88       include 'COMMON.THREAD'
  89       include 'COMMON.CONTACTS'
  90       include 'COMMON.TIME1'
  91 #ifdef MPL
  92       include 'COMMON.INFO'
  93 #endif
  94       character*4 sequence(maxres)
  95       integer rescode
  96       double precision x(maxvar)
  97       character*64 pdbfile
  98       dimension itype_pdb(maxres)
  99       logical seq_comp
 100 C
 101 C Body
 102 C
 103       if (indpdb.gt.0 .or. pdbref) then
 104         read(inp,'(a)') pdbfile
 105         write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 106         open(ipdbin,file=pdbfile,status='old',err=33)
 107         goto 34
 108   33    write (iout,'(a)') 'Error opening PDB file.'
 109         stop
 110   34    continue
 111         print *,'Begin reading pdb data'
 112         call readpdb
 113         print *,'Finished reading pdb data'
 114         write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 115         do i=1,nres
 116           itype_pdb(i)=itype(i)
 117         enddo
 118         close (ipdbin)
 119       endif
 120       if (indpdb.eq.0) then
 121 C Read sequence if not taken from the pdb file.
 122         if (iscode.gt.0) then
 123           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 124         else
 125           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 126         endif
 127 C Convert sequence to numeric code
 128         do i=1,nres
 129           itype(i)=rescode(i,sequence(i),iscode)
 130         enddo
 131       endif
 132       print *,nres
 133       print '(20i4)',(itype(i),i=1,nres)
 134       do i=1,nres
 135         if (itype(i).ne.20) then
 136           itel(i)=1
 137         else
 138           itel(i)=2
 139         endif
 140       enddo
 141       print *,'Call Read_Bridge.'
 142       call read_bridge
 143       nnt=1
 144       nct=nres
 145       print *,'NNT=',NNT,' NCT=',NCT
 146       if (itype(1).eq.21) nnt=2
 147       if (itype(nres).eq.21) nct=nct-1
 148       if (pdbref) then
 149         write (iout,'(a,i3)') 'nsup=',nsup
 150         nstart_seq=nnt
 151         if (nsup.le.(nct-nnt+1)) then
 152           do i=0,nct-nnt+1-nsup
 153             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 154               nstart_seq=nnt+i
 155               goto 111
 156             endif
 157           enddo
 158           write (iout,'(a)')
 159      &            'Error - sequences to be superposed do not match.'
 160           stop
 161         else
 162           do i=0,nsup-(nct-nnt+1)
 163             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 164      &      then
 165               nstart_sup=nstart_sup+i
 166               nsup=nct-nnt+1
 167               goto 111
 168             endif
 169           enddo
 170           write (iout,'(a)')
 171      &            'Error - sequences to be superposed do not match.'
 172         endif
 173   111   continue
 174         write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 175      &                 ' nstart_seq=',nstart_seq
 176       endif
 177       if (refstr) then
 178         if (.not.pdbref) then
 179           call read_angles(inp,*38)
 180           goto 39
 181    38     write (iout,'(a)') 'Error reading reference structure.'
 182 #ifdef MPL
 183           call mp_stopall(Error_Msg)
 184 #else
 185           stop 'Error reading reference structure'
 186 #endif
 187    39     call chainbuild
 188           nstart_sup=nnt
 189           nstart_seq=nnt
 190           nsup=nct-nnt+1
 191           do i=1,2*nres
 192             do j=1,3
 193               cref(j,i)=c(j,i)
 194             enddo
 195           enddo
 196         endif
 197         call contact(.true.,ncont_ref,icont_ref)
 198       endif
 199       return
 200       end
 201 c-----------------------------------------------------------------------------
 202       logical function seq_comp(itypea,itypeb,length)
 203       implicit none
 204       integer length,itypea(length),itypeb(length)
 205       integer i
 206       do i=1,length
 207         if (itypea(i).ne.itypeb(i)) then
 208           seq_comp=.false.
 209           return
 210         endif
 211       enddo
 212       seq_comp=.true.
 213       return
 214       end
 215 c-----------------------------------------------------------------------------
 216       subroutine read_bridge
 217 C Read information about disulfide bridges.
 218       include 'DIMENSIONS'
 219       include 'COMMON.IOUNITS'
 220       include 'COMMON.GEO'
 221       include 'COMMON.VAR'
 222       include 'COMMON.INTERACT'
 223       include 'COMMON.LOCAL'
 224       include 'COMMON.NAMES'
 225       include 'COMMON.CHAIN'
 226       include 'COMMON.FFIELD'
 227       include 'COMMON.SBRIDGE'
 228       include 'COMMON.HEADER'
 229       include 'COMMON.CONTROL'
 230       include 'COMMON.DBASE'
 231       include 'COMMON.THREAD'
 232       include 'COMMON.TIME1'
 233 #ifdef MPL
 234       include 'COMMON.INFO'
 235 #endif
 236 C Read bridging residues.
 237       read (inp,*) ns,(iss(i),i=1,ns)
 238       print *,'ns=',ns
 239 C Check whether the specified bridging residues are cystines.
 240       do i=1,ns
 241         if (itype(iss(i)).ne.1) then
 242           write (iout,'(2a,i3,a)')
 243      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 244      &   ' can form a disulfide bridge?!!!'
 245           write (*,'(2a,i3,a)')
 246      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 247      &   ' can form a disulfide bridge?!!!'
 248 #ifdef MPL
 249          call mp_stopall(error_msg)
 250 #else
 251          stop
 252 #endif
 253         endif
 254       enddo
 255 C Read preformed bridges.
 256       if (ns.gt.0) then
 257       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 258       if (nss.gt.0) then
 259         nhpb=nss
 260 C Check if the residues involved in bridges are in the specified list of
 261 C bridging residues.
 262         do i=1,nss
 263           do j=1,i-1
 264             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 265      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 266               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 267      &      ' contains residues present in other pairs.'
 268               write (*,'(a,i3,a)') 'Disulfide pair',i,
 269      &      ' contains residues present in other pairs.'
 270 #ifdef MPL
 271               call mp_stopall(error_msg)
 272 #else
 273               stop
 274 #endif
 275             endif
 276           enddo
 277           do j=1,ns
 278             if (ihpb(i).eq.iss(j)) goto 10
 279           enddo
 280           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 281    10     continue
 282           do j=1,ns
 283             if (jhpb(i).eq.iss(j)) goto 20
 284           enddo
 285           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 286    20     continue
 287           dhpb(i)=dbr
 288           forcon(i)=fbr
 289         enddo
 290         do i=1,nss
 291           ihpb(i)=ihpb(i)+nres
 292           jhpb(i)=jhpb(i)+nres
 293         enddo
 294       endif
 295       endif
 296       return
 297       end
 298 c----------------------------------------------------------------------------
 299       subroutine read_angles(kanal,*)
 300       include 'DIMENSIONS'
 301       include 'COMMON.GEO'
 302       include 'COMMON.VAR'
 303       include 'COMMON.CHAIN'
 304       include 'COMMON.IOUNITS'
 305       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 306       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 307       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 308       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 309       do i=1,nres
 310         theta(i)=deg2rad*theta(i)
 311         phi(i)=deg2rad*phi(i)
 312         alph(i)=deg2rad*alph(i)
 313         omeg(i)=deg2rad*omeg(i)
 314       enddo
 315       return
 316    10 return1
 317       end
 318 c----------------------------------------------------------------------------
 319       subroutine reada(rekord,lancuch,wartosc,default)
 320       implicit none
 321       character*(*) rekord,lancuch
 322       double precision wartosc,default
 323       integer ilen,iread
 324       external ilen
 325       iread=index(rekord,lancuch)
 326       if (iread.eq.0) then
 327         wartosc=default
 328         return
 329       endif
 330       iread=iread+ilen(lancuch)+1
 331       read (rekord(iread:),*) wartosc
 332       return
 333       end
 334 c----------------------------------------------------------------------------
 335       subroutine multreada(rekord,lancuch,tablica,dim,default)
 336       implicit none
 337       integer dim,i
 338       double precision tablica(dim),default
 339       character*(*) rekord,lancuch
 340       integer ilen,iread
 341       external ilen
 342       do i=1,dim
 343         tablica(i)=default
 344       enddo
 345       iread=index(rekord,lancuch)
 346       if (iread.eq.0) return
 347       iread=iread+ilen(lancuch)+1
 348       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 349    10 return
 350       end
 351 c----------------------------------------------------------------------------
 352       subroutine readi(rekord,lancuch,wartosc,default)
 353       implicit none
 354       character*(*) rekord,lancuch
 355       integer wartosc,default
 356       integer ilen,iread
 357       external ilen
 358       iread=index(rekord,lancuch)
 359       if (iread.eq.0) then
 360         wartosc=default
 361         return
 362       endif
 363       iread=iread+ilen(lancuch)+1
 364       read (rekord(iread:),*) wartosc
 365       return
 366       end
 367 c----------------------------------------------------------------------------
 368       subroutine split_string(rekord,tablica,dim,nsub)
 369       implicit none
 370       integer dim,nsub,i,ii,ll,kk
 371       character*(*) tablica(dim)
 372       character*(*) rekord
 373       integer ilen
 374       external ilen
 375       do i=1,dim
 376         tablica(i)=" "
 377       enddo
 378       ii=1
 379       ll = ilen(rekord)
 380       nsub=0
 381       do i=1,dim
 382 C Find the start of term name
 383         kk = 0
 384         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
 385           ii = ii+1
 386         enddo
 387 C Parse the name into TABLICA(i) until blank found
 388         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
 389           kk = kk+1
 390           tablica(i)(kk:kk)=rekord(ii:ii)
 391           ii = ii+1
 392         enddo
 393         if (kk.gt.0) nsub=nsub+1
 394         if (ii.gt.ll) return
 395       enddo
 396       return
 397       end
 398 c----------------------------------------------------------------------------
 399       subroutine card_concat(card)
 400       include 'DIMENSIONS'
 401       include 'COMMON.IOUNITS'
 402       character*(*) card
 403       character*80 karta,ucase
 404       external ilen
 405       read (inp,'(a)') karta
 406       karta=ucase(karta)
 407       card=' '
 408       do while (karta(80:80).eq.'&')
 409         card=card(:ilen(card)+1)//karta(:79)
 410         read (inp,'(a)') karta
 411         karta=ucase(karta)
 412       enddo
 413       card=card(:ilen(card)+1)//karta
 414       return
 415       end
 416 c----------------------------------------------------------------------------
 417       subroutine openunits
 418       include 'DIMENSIONS'
 419       include 'COMMON.IOUNITS'
 420       include 'COMMON.CONTROL'
 421       integer lenpre,lenpot,ilen
 422       external ilen
 423       character*16 cfrom_pdb,cprint
 424       character*16 ucase
 425       call getenv('INPUT',prefix)
 426       call getenv('OUTPUT',prefout)
 427       call getenv('INTIN',prefintin)
 428       call getenv('PDB',cfrom_pdb)
 429       call getenv('PRINTCOOR',cprint)
 430       from_cart = index(ucase(cfrom_pdb),'CART').gt.0
 431       from_cx = index(ucase(cfrom_pdb),'CX').gt.0
 432       lprint_cart = index(ucase(cprint),'CART').gt.0
 433       lprint_int = index(ucase(cprint),'INT').gt.0
 434 c      if (.not.lprint_cart .and. .not.lprint_int) lprint_int=.true.
 435       lenpre=ilen(prefix)
 436       lenout=ilen(prefout)
 437       lenint=ilen(prefintin)
 438 C Get the names and open the input files
 439       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 440       outname=prefout(:lenout)//'_clust.out'
 441       if (from_cart) then
 442         intinname=prefintin(:lenint)//'.x'
 443       else if (from_cx) then
 444         write (iout,*) "cx files: ",prefintin(:ilen(prefintin))
 445         call split_string(prefintin,cxfiles(1),maxfiles,nfiles)
 446         write (iout,*) "nfiles",nfiles
 447         write (iout,*) "Split cxfiles"
 448         do i=1,nfiles
 449           cxfiles(i)=cxfiles(i)(:ilen(cxfiles(i)))//'.cx'
 450           write (iout,*) cxfiles(i)(:ilen(cxfiles(i)))
 451         enddo
 452       else
 453         intinname=prefintin(:lenint)//'.int'
 454       endif
 455       if (lprint_cart) then
 456         if (.not. from_cx) then
 457           intname=prefintin(:lenint)//'_clust'//'.x'
 458         else
 459           intname=cxfiles(1)(:ilen(cxfiles(1)))//'_clust'//'.x'
 460         endif
 461       else if (lprint_int) then
 462         if (.not. from_cx) then
 463           intname=prefintin(:lenint)//'_clust'//'.int'
 464         else
 465           intname=cxfiles(1)(:ilen(cxfiles(1)))//'_clust'//'.int'
 466         endif
 467       endif
 468       rmsname=prefintin(:lenint)//'.rms'
 469       statinname=prefintin(:lenint)//'.stat'
 470       print *,statinname
 471       statoutname=prefintin(:lenint)//'_clust.stat'
 472       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 473      &       status='unknown')
 474       print *,'unit',jplot,' opened'
 475       if (.not. from_cx) open (intin,file=intinname,status='old')
 476       print *,'unit',intin,' opened'
 477       open (jrms,file=rmsname,status='unknown')
 478       open (jstatin,file=statinname,status='unknown')
 479       open (jstatout,file=statoutname,status='unknown')
 480 #ifdef AIX
 481       open(iout,file=outname,status='unknown')
 482       print *,'unit',iout,' opened'
 483       open(igeom,file=intname,status='unknown',position='append')
 484       print *,'unit',igeom,' opened'
 485 #elif defined(G77)
 486       open(iout,file=outname,status='unknown')
 487       open(igeom,file=intname,status='unknown',access='append')
 488 #else
 489       open(iout,file=outname,status='unknown')
 490       open(igeom,file=intname,status='unknown',position='append')
 491 #endif
 492       return
 493       end