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