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