source/cluster/wham/src/readrtns.F

   1       subroutine read_control
   2 C
   3 C Read molecular data
   4 C
   5       implicit none
   6       include 'DIMENSIONS'
   7       include 'sizesclu.dat'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.TIME1'
  10       include 'COMMON.SBRIDGE'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.HEADER'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.FREE'
  17       character*320 controlcard,ucase
  18 #ifdef MPL
  19       include 'COMMON.INFO'
  20 #endif
  21       integer i
  22
  23       read (INP,'(a80)') titel
  24       call card_concat(controlcard)
  25
  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)
  56       do i=1,nT
  57         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
  58       enddo
  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       if (min_var) iopt=1
  63       return
  64       end
  65 c--------------------------------------------------------------------------
  66       subroutine molread
  67 C
  68 C Read molecular data.
  69 C
  70       implicit none
  71       include 'DIMENSIONS'
  72       include 'COMMON.IOUNITS'
  73       include 'COMMON.GEO'
  74       include 'COMMON.VAR'
  75       include 'COMMON.INTERACT'
  76       include 'COMMON.LOCAL'
  77       include 'COMMON.NAMES'
  78       include 'COMMON.CHAIN'
  79       include 'COMMON.FFIELD'
  80       include 'COMMON.SBRIDGE'
  81       include 'COMMON.HEADER'
  82       include 'COMMON.CONTROL'
  83       include 'COMMON.CONTACTS'
  84       include 'COMMON.TIME1'
  85 #ifdef MPL
  86       include 'COMMON.INFO'
  87 #endif
  88       character*4 sequence(maxres)
  89       character*800 weightcard
  90       integer rescode
  91       double precision x(maxvar)
  92       integer itype_pdb(maxres)
  93       logical seq_comp
  94       integer i,j
  95 C
  96 C Body
  97 C
  98 C Read weights of the subsequent energy terms.
  99       call card_concat(weightcard)
 100       call reada(weightcard,'WSC',wsc,1.0d0)
 101       call reada(weightcard,'WLONG',wsc,wsc)
 102       call reada(weightcard,'WSCP',wscp,1.0d0)
 103       call reada(weightcard,'WELEC',welec,1.0D0)
 104       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 105       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 106       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 107       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 108       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 109       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 110       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 111       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 112       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 113       call reada(weightcard,'WBOND',wbond,1.0D0)
 114       call reada(weightcard,'WTOR',wtor,1.0D0)
 115       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 116       call reada(weightcard,'WANG',wang,1.0D0)
 117       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 118       call reada(weightcard,'SCAL14',scal14,0.4D0)
 119       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 120       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 121       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 122       if (index(weightcard,'SOFT').gt.0) ipot=6
 123 C 12/1/95 Added weight for the multi-body term WCORR
 124       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 125       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 126       weights(1)=wsc
 127       weights(2)=wscp
 128       weights(3)=welec
 129       weights(4)=wcorr
 130       weights(5)=wcorr5
 131       weights(6)=wcorr6
 132       weights(7)=wel_loc
 133       weights(8)=wturn3
 134       weights(9)=wturn4
 135       weights(10)=wturn6
 136       weights(11)=wang
 137       weights(12)=wscloc
 138       weights(13)=wtor
 139       weights(14)=wtor_d
 140       weights(15)=wstrain
 141       weights(16)=wvdwpp
 142       weights(17)=wbond
 143       weights(18)=scal14
 144       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 145      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 146      &  wturn4,wturn6
 147    10 format (/'Energy-term weights (unscaled):'//
 148      & 'WSCC=   ',f10.6,' (SC-SC)'/
 149      & 'WSCP=   ',f10.6,' (SC-p)'/
 150      & 'WELEC=  ',f10.6,' (p-p electr)'/
 151      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 152      & 'WBOND=  ',f10.6,' (stretching)'/
 153      & 'WANG=   ',f10.6,' (bending)'/
 154      & 'WSCLOC= ',f10.6,' (SC local)'/
 155      & 'WTOR=   ',f10.6,' (torsional)'/
 156      & 'WTORD=  ',f10.6,' (double torsional)'/
 157      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 158      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 159      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 160      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 161      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 162      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 163      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 164      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 165       if (wcorr4.gt.0.0d0) then
 166         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 167      &   'between contact pairs of peptide groups'
 168         write (iout,'(2(a,f5.3/))')
 169      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 170      &  'Range of quenching the correlation terms:',2*delt_corr
 171       else if (wcorr.gt.0.0d0) then
 172         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 173      &   'between contact pairs of peptide groups'
 174       endif
 175       write (iout,'(a,f8.3)')
 176      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 177       write (iout,'(a,f8.3)')
 178      &  'General scaling factor of SC-p interactions:',scalscp
 179       r0_corr=cutoff_corr-delt_corr
 180       do i=1,20
 181         aad(i,1)=scalscp*aad(i,1)
 182         aad(i,2)=scalscp*aad(i,2)
 183         bad(i,1)=scalscp*bad(i,1)
 184         bad(i,2)=scalscp*bad(i,2)
 185       enddo
 186
 187       call flush(iout)
 188       print *,'indpdb=',indpdb,' pdbref=',pdbref
 189
 190 C Read sequence if not taken from the pdb file.
 191       if (iscode.gt.0) then
 192         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 193       else
 194         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 195       endif
 196 C Convert sequence to numeric code
 197       do i=1,nres
 198         itype(i)=rescode(i,sequence(i),iscode)
 199       enddo
 200       print *,nres
 201       print '(20i4)',(itype(i),i=1,nres)
 202
 203       do i=1,nres
 204 #ifdef PROCOR
 205         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 206 #else
 207         if (itype(i).eq.21) then
 208 #endif
 209           itel(i)=0
 210 #ifdef PROCOR
 211         else if (itype(i+1).ne.20) then
 212 #else
 213         else if (itype(i).ne.20) then
 214 #endif
 215           itel(i)=1
 216         else
 217           itel(i)=2
 218         endif
 219       enddo
 220       write (iout,*) "ITEL"
 221       do i=1,nres-1
 222         write (iout,*) i,itype(i),itel(i)
 223       enddo
 224
 225       print *,'Call Read_Bridge.'
 226       call read_bridge
 227       nnt=1
 228       nct=nres
 229       print *,'NNT=',NNT,' NCT=',NCT
 230       if (itype(1).eq.21) nnt=2
 231       if (itype(nres).eq.21) nct=nct-1
 232       if (nstart.lt.nnt) nstart=nnt
 233       if (nend.gt.nct .or. nend.eq.0) nend=nct
 234       write (iout,*) "nstart",nstart," nend",nend
 235       nres0=nres
 236 c      if (pdbref) then
 237 c        read(inp,'(a)') pdbfile
 238 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 239 c        open(ipdbin,file=pdbfile,status='old',err=33)
 240 c        goto 34
 241 c  33    write (iout,'(a)') 'Error opening PDB file.'
 242 c        stop
 243 c  34    continue
 244 c        print *,'Begin reading pdb data'
 245 c        call readpdb
 246 c        print *,'Finished reading pdb data'
 247 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 248 c        do i=1,nres
 249 c          itype_pdb(i)=itype(i)
 250 c        enddo
 251 c        close (ipdbin)
 252 c        write (iout,'(a,i3)') 'nsup=',nsup
 253 c        nstart_seq=nnt
 254 c        if (nsup.le.(nct-nnt+1)) then
 255 c          do i=0,nct-nnt+1-nsup
 256 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 257 c              nstart_seq=nnt+i
 258 c              goto 111
 259 c            endif
 260 c          enddo
 261 c          write (iout,'(a)')
 262 c     &            'Error - sequences to be superposed do not match.'
 263 c          stop
 264 c        else
 265 c          do i=0,nsup-(nct-nnt+1)
 266 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 267 c     &      then
 268 c              nstart_sup=nstart_sup+i
 269 c              nsup=nct-nnt+1
 270 c              goto 111
 271 c            endif
 272 c          enddo
 273 c          write (iout,'(a)')
 274 c     &            'Error - sequences to be superposed do not match.'
 275 c        endif
 276 c  111   continue
 277 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 278 c     &                 ' nstart_seq=',nstart_seq
 279 c      endif
 280       call init_int_table
 281       call setup_var
 282       write (iout,*) "molread: REFSTR",refstr
 283       if (refstr) then
 284         if (.not.pdbref) then
 285           call read_angles(inp,*38)
 286           goto 39
 287    38     write (iout,'(a)') 'Error reading reference structure.'
 288 #ifdef MPL
 289           call mp_stopall(Error_Msg)
 290 #else
 291           stop 'Error reading reference structure'
 292 #endif
 293    39     call chainbuild
 294           nstart_sup=nnt
 295           nstart_seq=nnt
 296           nsup=nct-nnt+1
 297           do i=1,2*nres
 298             do j=1,3
 299               cref(j,i)=c(j,i)
 300             enddo
 301           enddo
 302         endif
 303         call contact(.true.,ncont_ref,icont_ref)
 304       endif
 305       return
 306       end
 307 c-----------------------------------------------------------------------------
 308       logical function seq_comp(itypea,itypeb,length)
 309       implicit none
 310       integer length,itypea(length),itypeb(length)
 311       integer i
 312       do i=1,length
 313         if (itypea(i).ne.itypeb(i)) then
 314           seq_comp=.false.
 315           return
 316         endif
 317       enddo
 318       seq_comp=.true.
 319       return
 320       end
 321 c-----------------------------------------------------------------------------
 322       subroutine read_bridge
 323 C Read information about disulfide bridges.
 324       implicit none
 325       include 'DIMENSIONS'
 326       include 'COMMON.IOUNITS'
 327       include 'COMMON.GEO'
 328       include 'COMMON.VAR'
 329       include 'COMMON.INTERACT'
 330       include 'COMMON.LOCAL'
 331       include 'COMMON.NAMES'
 332       include 'COMMON.CHAIN'
 333       include 'COMMON.FFIELD'
 334       include 'COMMON.SBRIDGE'
 335       include 'COMMON.HEADER'
 336       include 'COMMON.CONTROL'
 337       include 'COMMON.TIME1'
 338 #ifdef MPL
 339       include 'COMMON.INFO'
 340 #endif
 341       integer i,j
 342 C Read bridging residues.
 343       read (inp,*) ns,(iss(i),i=1,ns)
 344       print *,'ns=',ns
 345 C Check whether the specified bridging residues are cystines.
 346       do i=1,ns
 347         if (itype(iss(i)).ne.1) then
 348           write (iout,'(2a,i3,a)')
 349      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 350      &   ' can form a disulfide bridge?!!!'
 351           write (*,'(2a,i3,a)')
 352      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 353      &   ' can form a disulfide bridge?!!!'
 354 #ifdef MPL
 355          call mp_stopall(error_msg)
 356 #else
 357          stop
 358 #endif
 359         endif
 360       enddo
 361 C Read preformed bridges.
 362       if (ns.gt.0) then
 363       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 364       if (nss.gt.0) then
 365         nhpb=nss
 366 C Check if the residues involved in bridges are in the specified list of
 367 C bridging residues.
 368         do i=1,nss
 369           do j=1,i-1
 370             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 371      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 372               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 373      &      ' contains residues present in other pairs.'
 374               write (*,'(a,i3,a)') 'Disulfide pair',i,
 375      &      ' contains residues present in other pairs.'
 376 #ifdef MPL
 377               call mp_stopall(error_msg)
 378 #else
 379               stop
 380 #endif
 381             endif
 382           enddo
 383           do j=1,ns
 384             if (ihpb(i).eq.iss(j)) goto 10
 385           enddo
 386           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 387    10     continue
 388           do j=1,ns
 389             if (jhpb(i).eq.iss(j)) goto 20
 390           enddo
 391           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 392    20     continue
 393           dhpb(i)=dbr
 394           forcon(i)=fbr
 395         enddo
 396         do i=1,nss
 397           ihpb(i)=ihpb(i)+nres
 398           jhpb(i)=jhpb(i)+nres
 399         enddo
 400       endif
 401       endif
 402       return
 403       end
 404 c----------------------------------------------------------------------------
 405       subroutine read_angles(kanal,*)
 406       implicit none
 407       include 'DIMENSIONS'
 408       include 'COMMON.GEO'
 409       include 'COMMON.VAR'
 410       include 'COMMON.CHAIN'
 411       include 'COMMON.IOUNITS'
 412       integer i,kanal
 413       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 414       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 415       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 416       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 417       do i=1,nres
 418         theta(i)=deg2rad*theta(i)
 419         phi(i)=deg2rad*phi(i)
 420         alph(i)=deg2rad*alph(i)
 421         omeg(i)=deg2rad*omeg(i)
 422       enddo
 423       return
 424    10 return1
 425       end
 426 c----------------------------------------------------------------------------
 427       subroutine reada(rekord,lancuch,wartosc,default)
 428       implicit none
 429       character*(*) rekord,lancuch
 430       double precision wartosc,default
 431       integer ilen,iread
 432       external ilen
 433       iread=index(rekord,lancuch)
 434       if (iread.eq.0) then
 435         wartosc=default
 436         return
 437       endif
 438       iread=iread+ilen(lancuch)+1
 439       read (rekord(iread:),*) wartosc
 440       return
 441       end
 442 c----------------------------------------------------------------------------
 443       subroutine multreada(rekord,lancuch,tablica,dim,default)
 444       implicit none
 445       integer dim,i
 446       double precision tablica(dim),default
 447       character*(*) rekord,lancuch
 448       integer ilen,iread
 449       external ilen
 450       do i=1,dim
 451         tablica(i)=default
 452       enddo
 453       iread=index(rekord,lancuch)
 454       if (iread.eq.0) return
 455       iread=iread+ilen(lancuch)+1
 456       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 457    10 return
 458       end
 459 c----------------------------------------------------------------------------
 460       subroutine readi(rekord,lancuch,wartosc,default)
 461       implicit none
 462       character*(*) rekord,lancuch
 463       integer wartosc,default
 464       integer ilen,iread
 465       external ilen
 466       iread=index(rekord,lancuch)
 467       if (iread.eq.0) then
 468         wartosc=default
 469         return
 470       endif
 471       iread=iread+ilen(lancuch)+1
 472       read (rekord(iread:),*) wartosc
 473       return
 474       end
 475 c----------------------------------------------------------------------------
 476       subroutine card_concat(card)
 477       include 'DIMENSIONS'
 478       include 'COMMON.IOUNITS'
 479       character*(*) card
 480       character*80 karta,ucase
 481       external ilen
 482       read (inp,'(a)') karta
 483       karta=ucase(karta)
 484       card=' '
 485       do while (karta(80:80).eq.'&')
 486         card=card(:ilen(card)+1)//karta(:79)
 487         read (inp,'(a)') karta
 488         karta=ucase(karta)
 489       enddo
 490       card=card(:ilen(card)+1)//karta
 491       return
 492       end
 493 c----------------------------------------------------------------------------
 494       subroutine openunits
 495       implicit none
 496       include 'DIMENSIONS'
 497 #ifdef MPI
 498       include "mpif.h"
 499       character*3 liczba
 500       include "COMMON.MPI"
 501 #endif
 502       include 'COMMON.IOUNITS'
 503       include 'COMMON.CONTROL'
 504       integer lenpre,lenpot,ilen
 505       external ilen
 506       character*16 cformat,cprint
 507       character*16 ucase
 508       integer lenint,lenout
 509       call getenv('INPUT',prefix)
 510       call getenv('OUTPUT',prefout)
 511       call getenv('INTIN',prefintin)
 512       call getenv('COORD',cformat)
 513       call getenv('PRINTCOOR',cprint)
 514       call getenv('SCRATCHDIR',scratchdir)
 515       from_bx=.true.
 516       from_cx=.false.
 517       if (index(ucase(cformat),'CX').gt.0) then
 518         from_cx=.true.
 519         from_bx=.false.
 520       endif
 521       from_cart=.true.
 522       lenpre=ilen(prefix)
 523       lenout=ilen(prefout)
 524       lenint=ilen(prefintin)
 525 C Get the names and open the input files
 526       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 527 #ifdef MPI
 528       write (liczba,'(bz,i3.3)') me
 529       outname=prefout(:lenout)//'_clust.out_'//liczba
 530 #else
 531       outname=prefout(:lenout)//'_clust.out'
 532 #endif
 533       if (from_bx) then
 534         intinname=prefintin(:lenint)//'.bx'
 535       else if (from_cx) then
 536         intinname=prefintin(:lenint)//'.cx'
 537       else
 538         intinname=prefintin(:lenint)//'.int'
 539       endif
 540       rmsname=prefintin(:lenint)//'.rms'
 541       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 542      &       status='unknown')
 543       open (jrms,file=rmsname,status='unknown')
 544       open(iout,file=outname,status='unknown')
 545 C Get parameter filenames and open the parameter files.
 546       call getenv('BONDPAR',bondname)
 547       open (ibond,file=bondname,status='old')
 548       call getenv('THETPAR',thetname)
 549       open (ithep,file=thetname,status='old')
 550       call getenv('ROTPAR',rotname)
 551       open (irotam,file=rotname,status='old')
 552       call getenv('TORPAR',torname)
 553       open (itorp,file=torname,status='old')
 554       call getenv('TORDPAR',tordname)
 555       open (itordp,file=tordname,status='old')
 556       call getenv('FOURIER',fouriername)
 557       open (ifourier,file=fouriername,status='old')
 558       call getenv('ELEPAR',elename)
 559       open (ielep,file=elename,status='old')
 560       call getenv('SIDEPAR',sidename)
 561       open (isidep,file=sidename,status='old')
 562       call getenv('SIDEP',sidepname)
 563       open (isidep1,file=sidepname,status="old")
 564       call getenv('SCCORPAR',sccorname)
 565       open (isccor,file=sccorname,status="old")
 566 #ifndef OLDSCP
 567 C
 568 C 8/9/01 In the newest version SCp interaction constants are read from a file
 569 C Use -DOLDSCP to use hard-coded constants instead.
 570 C
 571       call getenv('SCPPAR',scpname)
 572       open (iscpp,file=scpname,status='old')
 573 #endif
 574       return
 575       end