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