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       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  63       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  64       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  65      & " CONSTR_DIST",constr_dist
  66       call flush(iout)
  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       implicit none
  76       include 'DIMENSIONS'
  77       include 'COMMON.IOUNITS'
  78       include 'COMMON.GEO'
  79       include 'COMMON.VAR'
  80       include 'COMMON.INTERACT'
  81       include 'COMMON.LOCAL'
  82       include 'COMMON.NAMES'
  83       include 'COMMON.CHAIN'
  84       include 'COMMON.FFIELD'
  85       include 'COMMON.SBRIDGE'
  86       include 'COMMON.HEADER'
  87       include 'COMMON.CONTROL'
  88       include 'COMMON.CONTACTS'
  89       include 'COMMON.TIME1'
  90       include 'COMMON.TORCNSTR'
  91 #ifdef MPL
  92       include 'COMMON.INFO'
  93 #endif
  94       character*4 sequence(maxres)
  95       character*800 weightcard
  96       integer rescode
  97       double precision x(maxvar)
  98       integer itype_pdb(maxres)
  99       logical seq_comp
 100       integer i,j
 101 C
 102 C Body
 103 C
 104 C Read weights of the subsequent energy terms.
 105       call card_concat(weightcard)
 106       call reada(weightcard,'WSC',wsc,1.0d0)
 107       call reada(weightcard,'WLONG',wsc,wsc)
 108       call reada(weightcard,'WSCP',wscp,1.0d0)
 109       call reada(weightcard,'WELEC',welec,1.0D0)
 110       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 111       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 112       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 113       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 114       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 115       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 116       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 117       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 118       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 119       call reada(weightcard,'WBOND',wbond,1.0D0)
 120       call reada(weightcard,'WTOR',wtor,1.0D0)
 121       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 122       call reada(weightcard,'WANG',wang,1.0D0)
 123       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 124       call reada(weightcard,'SCAL14',scal14,0.4D0)
 125       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 126       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 127       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 128       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 129       if (index(weightcard,'SOFT').gt.0) ipot=6
 130 C 12/1/95 Added weight for the multi-body term WCORR
 131       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 132       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 133       weights(1)=wsc
 134       weights(2)=wscp
 135       weights(3)=welec
 136       weights(4)=wcorr
 137       weights(5)=wcorr5
 138       weights(6)=wcorr6
 139       weights(7)=wel_loc
 140       weights(8)=wturn3
 141       weights(9)=wturn4
 142       weights(10)=wturn6
 143       weights(11)=wang
 144       weights(12)=wscloc
 145       weights(13)=wtor
 146       weights(14)=wtor_d
 147       weights(15)=wstrain
 148       weights(16)=wvdwpp
 149       weights(17)=wbond
 150       weights(18)=scal14
 151       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 152      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 153      &  wturn4,wturn6,wsccor
 154    10 format (/'Energy-term weights (unscaled):'//
 155      & 'WSCC=   ',f10.6,' (SC-SC)'/
 156      & 'WSCP=   ',f10.6,' (SC-p)'/
 157      & 'WELEC=  ',f10.6,' (p-p electr)'/
 158      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 159      & 'WBOND=  ',f10.6,' (stretching)'/
 160      & 'WANG=   ',f10.6,' (bending)'/
 161      & 'WSCLOC= ',f10.6,' (SC local)'/
 162      & 'WTOR=   ',f10.6,' (torsional)'/
 163      & 'WTORD=  ',f10.6,' (double torsional)'/
 164      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 165      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 166      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 167      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 168      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 169      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 170      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 171      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 172      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 173       if (wcorr4.gt.0.0d0) then
 174         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 175      &   'between contact pairs of peptide groups'
 176         write (iout,'(2(a,f5.3/))')
 177      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 178      &  'Range of quenching the correlation terms:',2*delt_corr
 179       else if (wcorr.gt.0.0d0) then
 180         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 181      &   'between contact pairs of peptide groups'
 182       endif
 183       write (iout,'(a,f8.3)')
 184      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 185       write (iout,'(a,f8.3)')
 186      &  'General scaling factor of SC-p interactions:',scalscp
 187       r0_corr=cutoff_corr-delt_corr
 188       do i=1,20
 189         aad(i,1)=scalscp*aad(i,1)
 190         aad(i,2)=scalscp*aad(i,2)
 191         bad(i,1)=scalscp*bad(i,1)
 192         bad(i,2)=scalscp*bad(i,2)
 193       enddo
 194
 195       call flush(iout)
 196       print *,'indpdb=',indpdb,' pdbref=',pdbref
 197
 198 C Read sequence if not taken from the pdb file.
 199       if (iscode.gt.0) then
 200         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 201       else
 202         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 203       endif
 204 C Convert sequence to numeric code
 205       do i=1,nres
 206         itype(i)=rescode(i,sequence(i),iscode)
 207       enddo
 208       print *,nres
 209       print '(20i4)',(itype(i),i=1,nres)
 210
 211       do i=1,nres
 212 #ifdef PROCOR
 213         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 214 #else
 215         if (itype(i).eq.21) then
 216 #endif
 217           itel(i)=0
 218 #ifdef PROCOR
 219         else if (itype(i+1).ne.20) then
 220 #else
 221         else if (itype(i).ne.20) then
 222 #endif
 223           itel(i)=1
 224         else
 225           itel(i)=2
 226         endif
 227       enddo
 228       write (iout,*) "ITEL"
 229       do i=1,nres-1
 230         write (iout,*) i,itype(i),itel(i)
 231       enddo
 232
 233       print *,'Call Read_Bridge.'
 234       call read_bridge
 235       if (with_dihed_constr) then
 236
 237       read (inp,*) ndih_constr
 238       if (ndih_constr.gt.0) then
 239         read (inp,*) ftors
 240         write (iout,*) 'FTORS',ftors
 241         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 242         write (iout,*)
 243      &   'There are',ndih_constr,' constraints on phi angles.'
 244         do i=1,ndih_constr
 245           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 246         enddo
 247         do i=1,ndih_constr
 248           phi0(i)=deg2rad*phi0(i)
 249           drange(i)=deg2rad*drange(i)
 250         enddo
 251       endif
 252
 253       endif
 254       nnt=1
 255       nct=nres
 256       print *,'NNT=',NNT,' NCT=',NCT
 257       if (itype(1).eq.21) nnt=2
 258       if (itype(nres).eq.21) nct=nct-1
 259       if (nstart.lt.nnt) nstart=nnt
 260       if (nend.gt.nct .or. nend.eq.0) nend=nct
 261       write (iout,*) "nstart",nstart," nend",nend
 262       nres0=nres
 263 c      if (pdbref) then
 264 c        read(inp,'(a)') pdbfile
 265 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 266 c        open(ipdbin,file=pdbfile,status='old',err=33)
 267 c        goto 34
 268 c  33    write (iout,'(a)') 'Error opening PDB file.'
 269 c        stop
 270 c  34    continue
 271 c        print *,'Begin reading pdb data'
 272 c        call readpdb
 273 c        print *,'Finished reading pdb data'
 274 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 275 c        do i=1,nres
 276 c          itype_pdb(i)=itype(i)
 277 c        enddo
 278 c        close (ipdbin)
 279 c        write (iout,'(a,i3)') 'nsup=',nsup
 280 c        nstart_seq=nnt
 281 c        if (nsup.le.(nct-nnt+1)) then
 282 c          do i=0,nct-nnt+1-nsup
 283 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 284 c              nstart_seq=nnt+i
 285 c              goto 111
 286 c            endif
 287 c          enddo
 288 c          write (iout,'(a)')
 289 c     &            'Error - sequences to be superposed do not match.'
 290 c          stop
 291 c        else
 292 c          do i=0,nsup-(nct-nnt+1)
 293 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 294 c     &      then
 295 c              nstart_sup=nstart_sup+i
 296 c              nsup=nct-nnt+1
 297 c              goto 111
 298 c            endif
 299 c          enddo
 300 c          write (iout,'(a)')
 301 c     &            'Error - sequences to be superposed do not match.'
 302 c        endif
 303 c  111   continue
 304 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 305 c     &                 ' nstart_seq=',nstart_seq
 306 c      endif
 307       call init_int_table
 308       call setup_var
 309       write (iout,*) "molread: REFSTR",refstr
 310       if (refstr) then
 311         if (.not.pdbref) then
 312           call read_angles(inp,*38)
 313           goto 39
 314    38     write (iout,'(a)') 'Error reading reference structure.'
 315 #ifdef MPL
 316           call mp_stopall(Error_Msg)
 317 #else
 318           stop 'Error reading reference structure'
 319 #endif
 320    39     call chainbuild
 321           nstart_sup=nnt
 322           nstart_seq=nnt
 323           nsup=nct-nnt+1
 324           do i=1,2*nres
 325             do j=1,3
 326               cref(j,i)=c(j,i)
 327             enddo
 328           enddo
 329         endif
 330         call contact(.true.,ncont_ref,icont_ref)
 331       endif
 332 c Read distance restraints
 333       if (constr_dist.gt.0) then
 334         call read_dist_constr
 335         call hpb_partition
 336       endif
 337       return
 338       end
 339 c-----------------------------------------------------------------------------
 340       logical function seq_comp(itypea,itypeb,length)
 341       implicit none
 342       integer length,itypea(length),itypeb(length)
 343       integer i
 344       do i=1,length
 345         if (itypea(i).ne.itypeb(i)) then
 346           seq_comp=.false.
 347           return
 348         endif
 349       enddo
 350       seq_comp=.true.
 351       return
 352       end
 353 c-----------------------------------------------------------------------------
 354       subroutine read_bridge
 355 C Read information about disulfide bridges.
 356       implicit none
 357       include 'DIMENSIONS'
 358       include 'COMMON.IOUNITS'
 359       include 'COMMON.GEO'
 360       include 'COMMON.VAR'
 361       include 'COMMON.INTERACT'
 362       include 'COMMON.LOCAL'
 363       include 'COMMON.NAMES'
 364       include 'COMMON.CHAIN'
 365       include 'COMMON.FFIELD'
 366       include 'COMMON.SBRIDGE'
 367       include 'COMMON.HEADER'
 368       include 'COMMON.CONTROL'
 369       include 'COMMON.TIME1'
 370 #ifdef MPL
 371       include 'COMMON.INFO'
 372 #endif
 373       integer i,j
 374 C Read bridging residues.
 375       read (inp,*) ns,(iss(i),i=1,ns)
 376       print *,'ns=',ns
 377 C Check whether the specified bridging residues are cystines.
 378       do i=1,ns
 379         if (itype(iss(i)).ne.1) then
 380           write (iout,'(2a,i3,a)')
 381      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 382      &   ' can form a disulfide bridge?!!!'
 383           write (*,'(2a,i3,a)')
 384      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 385      &   ' can form a disulfide bridge?!!!'
 386 #ifdef MPL
 387          call mp_stopall(error_msg)
 388 #else
 389          stop
 390 #endif
 391         endif
 392       enddo
 393 C Read preformed bridges.
 394       if (ns.gt.0) then
 395       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 396       if (nss.gt.0) then
 397         nhpb=nss
 398 C Check if the residues involved in bridges are in the specified list of
 399 C bridging residues.
 400         do i=1,nss
 401           do j=1,i-1
 402             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 403      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 404               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 405      &      ' contains residues present in other pairs.'
 406               write (*,'(a,i3,a)') 'Disulfide pair',i,
 407      &      ' contains residues present in other pairs.'
 408 #ifdef MPL
 409               call mp_stopall(error_msg)
 410 #else
 411               stop
 412 #endif
 413             endif
 414           enddo
 415           do j=1,ns
 416             if (ihpb(i).eq.iss(j)) goto 10
 417           enddo
 418           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 419    10     continue
 420           do j=1,ns
 421             if (jhpb(i).eq.iss(j)) goto 20
 422           enddo
 423           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 424    20     continue
 425 c          dhpb(i)=dbr
 426 c          forcon(i)=fbr
 427         enddo
 428         do i=1,nss
 429           ihpb(i)=ihpb(i)+nres
 430           jhpb(i)=jhpb(i)+nres
 431         enddo
 432       endif
 433       endif
 434       return
 435       end
 436 c----------------------------------------------------------------------------
 437       subroutine read_angles(kanal,*)
 438       implicit none
 439       include 'DIMENSIONS'
 440       include 'COMMON.GEO'
 441       include 'COMMON.VAR'
 442       include 'COMMON.CHAIN'
 443       include 'COMMON.IOUNITS'
 444       integer i,kanal
 445       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 446       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 447       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 448       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 449       do i=1,nres
 450         theta(i)=deg2rad*theta(i)
 451         phi(i)=deg2rad*phi(i)
 452         alph(i)=deg2rad*alph(i)
 453         omeg(i)=deg2rad*omeg(i)
 454       enddo
 455       return
 456    10 return1
 457       end
 458 c----------------------------------------------------------------------------
 459       subroutine reada(rekord,lancuch,wartosc,default)
 460       implicit none
 461       character*(*) rekord,lancuch
 462       double precision wartosc,default
 463       integer ilen,iread
 464       external ilen
 465       iread=index(rekord,lancuch)
 466       if (iread.eq.0) then
 467         wartosc=default
 468         return
 469       endif
 470       iread=iread+ilen(lancuch)+1
 471       read (rekord(iread:),*) wartosc
 472       return
 473       end
 474 c----------------------------------------------------------------------------
 475       subroutine multreada(rekord,lancuch,tablica,dim,default)
 476       implicit none
 477       integer dim,i
 478       double precision tablica(dim),default
 479       character*(*) rekord,lancuch
 480       integer ilen,iread
 481       external ilen
 482       do i=1,dim
 483         tablica(i)=default
 484       enddo
 485       iread=index(rekord,lancuch)
 486       if (iread.eq.0) return
 487       iread=iread+ilen(lancuch)+1
 488       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 489    10 return
 490       end
 491 c----------------------------------------------------------------------------
 492       subroutine readi(rekord,lancuch,wartosc,default)
 493       implicit none
 494       character*(*) rekord,lancuch
 495       integer wartosc,default
 496       integer ilen,iread
 497       external ilen
 498       iread=index(rekord,lancuch)
 499       if (iread.eq.0) then
 500         wartosc=default
 501         return
 502       endif
 503       iread=iread+ilen(lancuch)+1
 504       read (rekord(iread:),*) wartosc
 505       return
 506       end
 507 c----------------------------------------------------------------------------
 508       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 509       implicit none
 510       integer dim,i
 511       integer tablica(dim),default
 512       character*(*) rekord,lancuch
 513       character*80 aux
 514       integer ilen,iread
 515       external ilen
 516       do i=1,dim
 517         tablica(i)=default
 518       enddo
 519       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 520       if (iread.eq.0) return
 521       iread=iread+ilen(lancuch)+1
 522       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 523    10 return
 524       end
 525 c----------------------------------------------------------------------------
 526       subroutine card_concat(card)
 527       include 'DIMENSIONS'
 528       include 'COMMON.IOUNITS'
 529       character*(*) card
 530       character*80 karta,ucase
 531       external ilen
 532       read (inp,'(a)') karta
 533       karta=ucase(karta)
 534       card=' '
 535       do while (karta(80:80).eq.'&')
 536         card=card(:ilen(card)+1)//karta(:79)
 537         read (inp,'(a)') karta
 538         karta=ucase(karta)
 539       enddo
 540       card=card(:ilen(card)+1)//karta
 541       return
 542       end
 543 c----------------------------------------------------------------------------
 544       subroutine openunits
 545       implicit none
 546       include 'DIMENSIONS'
 547 #ifdef MPI
 548       include "mpif.h"
 549       character*3 liczba
 550       include "COMMON.MPI"
 551 #endif
 552       include 'COMMON.IOUNITS'
 553       include 'COMMON.CONTROL'
 554       integer lenpre,lenpot,ilen
 555       external ilen
 556       character*16 cformat,cprint
 557       character*16 ucase
 558       integer lenint,lenout
 559       call getenv('INPUT',prefix)
 560       call getenv('OUTPUT',prefout)
 561       call getenv('INTIN',prefintin)
 562       call getenv('COORD',cformat)
 563       call getenv('PRINTCOOR',cprint)
 564       call getenv('SCRATCHDIR',scratchdir)
 565       from_bx=.true.
 566       from_cx=.false.
 567       if (index(ucase(cformat),'CX').gt.0) then
 568         from_cx=.true.
 569         from_bx=.false.
 570       endif
 571       from_cart=.true.
 572       lenpre=ilen(prefix)
 573       lenout=ilen(prefout)
 574       lenint=ilen(prefintin)
 575 C Get the names and open the input files
 576       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 577 #ifdef MPI
 578       write (liczba,'(bz,i3.3)') me
 579       outname=prefout(:lenout)//'_clust.out_'//liczba
 580 #else
 581       outname=prefout(:lenout)//'_clust.out'
 582 #endif
 583       if (from_bx) then
 584         intinname=prefintin(:lenint)//'.bx'
 585       else if (from_cx) then
 586         intinname=prefintin(:lenint)//'.cx'
 587       else
 588         intinname=prefintin(:lenint)//'.int'
 589       endif
 590       rmsname=prefintin(:lenint)//'.rms'
 591       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 592      &       status='unknown')
 593       open (jrms,file=rmsname,status='unknown')
 594       open(iout,file=outname,status='unknown')
 595 C Get parameter filenames and open the parameter files.
 596       call getenv('BONDPAR',bondname)
 597       open (ibond,file=bondname,status='old')
 598       call getenv('THETPAR',thetname)
 599       open (ithep,file=thetname,status='old')
 600       call getenv('ROTPAR',rotname)
 601       open (irotam,file=rotname,status='old')
 602       call getenv('TORPAR',torname)
 603       open (itorp,file=torname,status='old')
 604       call getenv('TORDPAR',tordname)
 605       open (itordp,file=tordname,status='old')
 606       call getenv('FOURIER',fouriername)
 607       open (ifourier,file=fouriername,status='old')
 608       call getenv('ELEPAR',elename)
 609       open (ielep,file=elename,status='old')
 610       call getenv('SIDEPAR',sidename)
 611       open (isidep,file=sidename,status='old')
 612       call getenv('SIDEP',sidepname)
 613       open (isidep1,file=sidepname,status="old")
 614       call getenv('SCCORPAR',sccorname)
 615       open (isccor,file=sccorname,status="old")
 616 #ifndef OLDSCP
 617 C
 618 C 8/9/01 In the newest version SCp interaction constants are read from a file
 619 C Use -DOLDSCP to use hard-coded constants instead.
 620 C
 621       call getenv('SCPPAR',scpname)
 622       open (iscpp,file=scpname,status='old')
 623 #endif
 624       return
 625       end
 626 c-------------------------------------------------------------------------------
 627       subroutine read_dist_constr
 628       implicit real*8 (a-h,o-z)
 629       include 'DIMENSIONS'
 630       include 'COMMON.CONTROL'
 631       include 'COMMON.CHAIN'
 632       include 'COMMON.IOUNITS'
 633       include 'COMMON.SBRIDGE'
 634       integer ifrag_(2,100),ipair_(2,100)
 635       double precision wfrag_(100),wpair_(100)
 636       character*500 controlcard
 637 c      write (iout,*) "Calling read_dist_constr"
 638 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 639       call card_concat(controlcard)
 640 c      call flush(iout)
 641 c      write (iout,'(a)') controlcard
 642       call readi(controlcard,"NFRAG",nfrag_,0)
 643       call readi(controlcard,"NPAIR",npair_,0)
 644       call readi(controlcard,"NDIST",ndist_,0)
 645       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 646       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 647       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 648       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 649       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 650       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 651       write (iout,*) "IFRAG"
 652       do i=1,nfrag_
 653         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 654       enddo
 655       write (iout,*) "IPAIR"
 656       do i=1,npair_
 657         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 658       enddo
 659       call flush(iout)
 660       if (.not.refstr .and. nfrag_.gt.0) then
 661         write (iout,*)
 662      &  "ERROR: no reference structure to compute distance restraints"
 663         write (iout,*)
 664      &  "Restraints must be specified explicitly (NDIST=number)"
 665         stop
 666       endif
 667       if (nfrag_.lt.2 .and. npair_.gt.0) then
 668         write (iout,*) "ERROR: Less than 2 fragments specified",
 669      &   " but distance restraints between pairs requested"
 670         stop
 671       endif
 672       call flush(iout)
 673       do i=1,nfrag_
 674         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 675         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 676      &    ifrag_(2,i)=nstart_sup+nsup-1
 677 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 678         call flush(iout)
 679         if (wfrag_(i).gt.0.0d0) then
 680         do j=ifrag_(1,i),ifrag_(2,i)-1
 681           do k=j+1,ifrag_(2,i)
 682             write (iout,*) "j",j," k",k
 683             ddjk=dist(j,k)
 684             if (constr_dist.eq.1) then
 685             nhpb=nhpb+1
 686             ihpb(nhpb)=j
 687             jhpb(nhpb)=k
 688               dhpb(nhpb)=ddjk
 689             forcon(nhpb)=wfrag_(i)
 690             else if (constr_dist.eq.2) then
 691               if (ddjk.le.dist_cut) then
 692                 nhpb=nhpb+1
 693                 ihpb(nhpb)=j
 694                 jhpb(nhpb)=k
 695                 dhpb(nhpb)=ddjk
 696                 forcon(nhpb)=wfrag_(i)
 697               endif
 698             else
 699               nhpb=nhpb+1
 700               ihpb(nhpb)=j
 701               jhpb(nhpb)=k
 702               dhpb(nhpb)=ddjk
 703               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 704             endif
 705             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 706      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 707           enddo
 708         enddo
 709         endif
 710       enddo
 711       do i=1,npair_
 712         if (wpair_(i).gt.0.0d0) then
 713         ii = ipair_(1,i)
 714         jj = ipair_(2,i)
 715         if (ii.gt.jj) then
 716           itemp=ii
 717           ii=jj
 718           jj=itemp
 719         endif
 720         do j=ifrag_(1,ii),ifrag_(2,ii)
 721           do k=ifrag_(1,jj),ifrag_(2,jj)
 722             nhpb=nhpb+1
 723             ihpb(nhpb)=j
 724             jhpb(nhpb)=k
 725             forcon(nhpb)=wpair_(i)
 726             dhpb(nhpb)=dist(j,k)
 727             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 728      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 729           enddo
 730         enddo
 731         endif
 732       enddo
 733       do i=1,ndist_
 734         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 735      &     ibecarb(i),forcon(nhpb+1)
 736         if (forcon(nhpb+1).gt.0.0d0) then
 737           nhpb=nhpb+1
 738           if (ibecarb(i).gt.0) then
 739             ihpb(i)=ihpb(i)+nres
 740             jhpb(i)=jhpb(i)+nres
 741           endif
 742           if (dhpb(nhpb).eq.0.0d0)
 743      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 744         endif
 745       enddo
 746       do i=1,nhpb
 747           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 748      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 749       enddo
 750       call flush(iout)
 751       return
 752       end