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         read (inp,*) ftors
 293         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),i=1,ndih_constr)
 296         write (iout,*)
 297      &   'There are',ndih_constr,' constraints on phi angles.'
 298         do i=1,ndih_constr
 299           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 300         enddo
 301         do i=1,ndih_constr
 302           phi0(i)=deg2rad*phi0(i)
 303           drange(i)=deg2rad*drange(i)
 304         enddo
 305       endif ! endif ndif_constr.gt.0
 306       endif ! with_dihed_constr
 307       nnt=1
 308       nct=nres
 309       print *,'NNT=',NNT,' NCT=',NCT
 310       if (itype(1).eq.ntyp1) nnt=2
 311       if (itype(nres).eq.ntyp1) nct=nct-1
 312       if (nstart.lt.nnt) nstart=nnt
 313       if (nend.gt.nct .or. nend.eq.0) nend=nct
 314       write (iout,*) "nstart",nstart," nend",nend
 315       nres0=nres
 316 c      if (pdbref) then
 317 c        read(inp,'(a)') pdbfile
 318 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 319 c        open(ipdbin,file=pdbfile,status='old',err=33)
 320 c        goto 34
 321 c  33    write (iout,'(a)') 'Error opening PDB file.'
 322 c        stop
 323 c  34    continue
 324 c        print *,'Begin reading pdb data'
 325 c        call readpdb
 326 c        print *,'Finished reading pdb data'
 327 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 328 c        do i=1,nres
 329 c          itype_pdb(i)=itype(i)
 330 c        enddo
 331 c        close (ipdbin)
 332 c        write (iout,'(a,i3)') 'nsup=',nsup
 333 c        nstart_seq=nnt
 334 c        if (nsup.le.(nct-nnt+1)) then
 335 c          do i=0,nct-nnt+1-nsup
 336 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 337 c              nstart_seq=nnt+i
 338 c              goto 111
 339 c            endif
 340 c          enddo
 341 c          write (iout,'(a)')
 342 c     &            'Error - sequences to be superposed do not match.'
 343 c          stop
 344 c        else
 345 c          do i=0,nsup-(nct-nnt+1)
 346 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 347 c     &      then
 348 c              nstart_sup=nstart_sup+i
 349 c              nsup=nct-nnt+1
 350 c              goto 111
 351 c            endif
 352 c          enddo
 353 c          write (iout,'(a)')
 354 c     &            'Error - sequences to be superposed do not match.'
 355 c        endif
 356 c  111   continue
 357 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 358 c     &                 ' nstart_seq=',nstart_seq
 359 c      endif
 360       call init_int_table
 361       call setup_var
 362       write (iout,*) "molread: REFSTR",refstr
 363       if (refstr) then
 364         if (.not.pdbref) then
 365           call read_angles(inp,*38)
 366           goto 39
 367    38     write (iout,'(a)') 'Error reading reference structure.'
 368 #ifdef MPL
 369           call mp_stopall(Error_Msg)
 370 #else
 371           stop 'Error reading reference structure'
 372 #endif
 373    39     call chainbuild
 374           nstart_sup=nnt
 375           nstart_seq=nnt
 376           nsup=nct-nnt+1
 377           kkk=1
 378           do i=1,2*nres
 379             do j=1,3
 380               cref(j,i,kkk)=c(j,i)
 381             enddo
 382           enddo
 383         endif
 384         call contact(.true.,ncont_ref,icont_ref)
 385       endif
 386        if (ns.gt.0) then
 387 C        write (iout,'(/a,i3,a)')
 388 C     &  'The chain contains',ns,' disulfide-bridging cysteines.'
 389         write (iout,'(20i4)') (iss(i),i=1,ns)
 390        if (dyn_ss) then
 391           write(iout,*)"Running with dynamic disulfide-bond formation"
 392        else
 393         write (iout,'(/a/)') 'Pre-formed links are:'
 394         do i=1,nss
 395           i1=ihpb(i)-nres
 396           i2=jhpb(i)-nres
 397           it1=itype(i1)
 398           it2=itype(i2)
 399           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 400      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 401      &    ebr,forcon(i)
 402         enddo
 403         write (iout,'(a)')
 404        endif
 405       endif
 406       if (ns.gt.0.and.dyn_ss) then
 407           do i=nss+1,nhpb
 408             ihpb(i-nss)=ihpb(i)
 409             jhpb(i-nss)=jhpb(i)
 410             forcon(i-nss)=forcon(i)
 411             dhpb(i-nss)=dhpb(i)
 412           enddo
 413           nhpb=nhpb-nss
 414           nss=0
 415           call hpb_partition
 416           do i=1,ns
 417             dyn_ss_mask(iss(i))=.true.
 418           enddo
 419       endif
 420 c Read distance restraints
 421       if (constr_dist.gt.0) then
 422         call read_dist_constr
 423         call hpb_partition
 424       endif
 425       return
 426       end
 427 c-----------------------------------------------------------------------------
 428       logical function seq_comp(itypea,itypeb,length)
 429       implicit none
 430       integer length,itypea(length),itypeb(length)
 431       integer i
 432       do i=1,length
 433         if (itypea(i).ne.itypeb(i)) then
 434           seq_comp=.false.
 435           return
 436         endif
 437       enddo
 438       seq_comp=.true.
 439       return
 440       end
 441 c-----------------------------------------------------------------------------
 442       subroutine read_bridge
 443 C Read information about disulfide bridges.
 444       implicit none
 445       include 'DIMENSIONS'
 446       include 'COMMON.IOUNITS'
 447       include 'COMMON.GEO'
 448       include 'COMMON.VAR'
 449       include 'COMMON.INTERACT'
 450       include 'COMMON.LOCAL'
 451       include 'COMMON.NAMES'
 452       include 'COMMON.CHAIN'
 453       include 'COMMON.FFIELD'
 454       include 'COMMON.SBRIDGE'
 455       include 'COMMON.HEADER'
 456       include 'COMMON.CONTROL'
 457       include 'COMMON.TIME1'
 458 #ifdef MPL
 459       include 'COMMON.INFO'
 460 #endif
 461       integer i,j
 462 C Read bridging residues.
 463       read (inp,*) ns,(iss(i),i=1,ns)
 464       print *,'ns=',ns
 465 C Check whether the specified bridging residues are cystines.
 466       do i=1,ns
 467         if (itype(iss(i)).ne.1) then
 468           write (iout,'(2a,i3,a)')
 469      &   'Do you REALLY think that the residue ',
 470      &    restyp(itype(iss(i))),i,
 471      &   ' can form a disulfide bridge?!!!'
 472           write (*,'(2a,i3,a)')
 473      &   'Do you REALLY think that the residue ',
 474      &   restyp(itype(iss(i))),i,
 475      &   ' can form a disulfide bridge?!!!'
 476 #ifdef MPL
 477          call mp_stopall(error_msg)
 478 #else
 479          stop
 480 #endif
 481         endif
 482       enddo
 483 C Read preformed bridges.
 484       if (ns.gt.0) then
 485       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 486       if (nss.gt.0) then
 487         nhpb=nss
 488 C Check if the residues involved in bridges are in the specified list of
 489 C bridging residues.
 490         do i=1,nss
 491           do j=1,i-1
 492             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 493      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 494               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 495      &      ' contains residues present in other pairs.'
 496               write (*,'(a,i3,a)') 'Disulfide pair',i,
 497      &      ' contains residues present in other pairs.'
 498 #ifdef MPL
 499               call mp_stopall(error_msg)
 500 #else
 501               stop
 502 #endif
 503             endif
 504           enddo
 505           do j=1,ns
 506             if (ihpb(i).eq.iss(j)) goto 10
 507           enddo
 508           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 509    10     continue
 510           do j=1,ns
 511             if (jhpb(i).eq.iss(j)) goto 20
 512           enddo
 513           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 514    20     continue
 515 C          dhpb(i)=dbr
 516 C          forcon(i)=fbr
 517         enddo
 518         do i=1,nss
 519           ihpb(i)=ihpb(i)+nres
 520           jhpb(i)=jhpb(i)+nres
 521         enddo
 522       endif
 523       endif
 524       return
 525       end
 526 c----------------------------------------------------------------------------
 527       subroutine read_angles(kanal,*)
 528       implicit none
 529       include 'DIMENSIONS'
 530       include 'COMMON.GEO'
 531       include 'COMMON.VAR'
 532       include 'COMMON.CHAIN'
 533       include 'COMMON.IOUNITS'
 534       integer i,kanal
 535       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 536       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 537       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 538       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 539       do i=1,nres
 540         theta(i)=deg2rad*theta(i)
 541         phi(i)=deg2rad*phi(i)
 542         alph(i)=deg2rad*alph(i)
 543         omeg(i)=deg2rad*omeg(i)
 544       enddo
 545       return
 546    10 return1
 547       end
 548 c----------------------------------------------------------------------------
 549       subroutine reada(rekord,lancuch,wartosc,default)
 550       implicit none
 551       character*(*) rekord,lancuch
 552       double precision wartosc,default
 553       integer ilen,iread
 554       external ilen
 555       iread=index(rekord,lancuch)
 556       if (iread.eq.0) then
 557         wartosc=default
 558         return
 559       endif
 560       iread=iread+ilen(lancuch)+1
 561       read (rekord(iread:),*) wartosc
 562       return
 563       end
 564 c----------------------------------------------------------------------------
 565       subroutine multreada(rekord,lancuch,tablica,dim,default)
 566       implicit none
 567       integer dim,i
 568       double precision tablica(dim),default
 569       character*(*) rekord,lancuch
 570       integer ilen,iread
 571       external ilen
 572       do i=1,dim
 573         tablica(i)=default
 574       enddo
 575       iread=index(rekord,lancuch)
 576       if (iread.eq.0) return
 577       iread=iread+ilen(lancuch)+1
 578       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 579    10 return
 580       end
 581 c----------------------------------------------------------------------------
 582       subroutine readi(rekord,lancuch,wartosc,default)
 583       implicit none
 584       character*(*) rekord,lancuch
 585       integer wartosc,default
 586       integer ilen,iread
 587       external ilen
 588       iread=index(rekord,lancuch)
 589       if (iread.eq.0) then
 590         wartosc=default
 591         return
 592       endif
 593       iread=iread+ilen(lancuch)+1
 594       read (rekord(iread:),*) wartosc
 595       return
 596       end
 597 C----------------------------------------------------------------------
 598       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 599       implicit none
 600       integer dim,i
 601       integer tablica(dim),default
 602       character*(*) rekord,lancuch
 603       character*80 aux
 604       integer ilen,iread
 605       external ilen
 606       do i=1,dim
 607         tablica(i)=default
 608       enddo
 609       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 610       if (iread.eq.0) return
 611       iread=iread+ilen(lancuch)+1
 612       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 613    10 return
 614       end
 615
 616 c----------------------------------------------------------------------------
 617       subroutine card_concat(card)
 618       include 'DIMENSIONS'
 619       include 'COMMON.IOUNITS'
 620       character*(*) card
 621       character*80 karta,ucase
 622       external ilen
 623       read (inp,'(a)') karta
 624       karta=ucase(karta)
 625       card=' '
 626       do while (karta(80:80).eq.'&')
 627         card=card(:ilen(card)+1)//karta(:79)
 628         read (inp,'(a)') karta
 629         karta=ucase(karta)
 630       enddo
 631       card=card(:ilen(card)+1)//karta
 632       return
 633       end
 634 c----------------------------------------------------------------------------
 635       subroutine openunits
 636       implicit none
 637       include 'DIMENSIONS'
 638 #ifdef MPI
 639       include "mpif.h"
 640       character*3 liczba
 641       include "COMMON.MPI"
 642 #endif
 643       include 'COMMON.IOUNITS'
 644       include 'COMMON.CONTROL'
 645       integer lenpre,lenpot,ilen
 646       external ilen
 647       character*16 cformat,cprint
 648       character*16 ucase
 649       integer lenint,lenout
 650       call getenv('INPUT',prefix)
 651       call getenv('OUTPUT',prefout)
 652       call getenv('INTIN',prefintin)
 653       call getenv('COORD',cformat)
 654       call getenv('PRINTCOOR',cprint)
 655       call getenv('SCRATCHDIR',scratchdir)
 656       from_bx=.true.
 657       from_cx=.false.
 658       if (index(ucase(cformat),'CX').gt.0) then
 659         from_cx=.true.
 660         from_bx=.false.
 661       endif
 662       from_cart=.true.
 663       lenpre=ilen(prefix)
 664       lenout=ilen(prefout)
 665       lenint=ilen(prefintin)
 666 C Get the names and open the input files
 667       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 668 #ifdef MPI
 669       write (liczba,'(bz,i3.3)') me
 670       outname=prefout(:lenout)//'_clust.out_'//liczba
 671 #else
 672       outname=prefout(:lenout)//'_clust.out'
 673 #endif
 674       if (from_bx) then
 675         intinname=prefintin(:lenint)//'.bx'
 676       else if (from_cx) then
 677         intinname=prefintin(:lenint)//'.cx'
 678       else
 679         intinname=prefintin(:lenint)//'.int'
 680       endif
 681       rmsname=prefintin(:lenint)//'.rms'
 682       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 683      &       status='unknown')
 684       open (jrms,file=rmsname,status='unknown')
 685       open(iout,file=outname,status='unknown')
 686 C Get parameter filenames and open the parameter files.
 687       call getenv('BONDPAR',bondname)
 688       open (ibond,file=bondname,status='old')
 689       call getenv('THETPAR',thetname)
 690       open (ithep,file=thetname,status='old')
 691       call getenv('ROTPAR',rotname)
 692       open (irotam,file=rotname,status='old')
 693       call getenv('TORPAR',torname)
 694       open (itorp,file=torname,status='old')
 695       call getenv('TORDPAR',tordname)
 696       open (itordp,file=tordname,status='old')
 697       call getenv('FOURIER',fouriername)
 698       open (ifourier,file=fouriername,status='old')
 699       call getenv('ELEPAR',elename)
 700       open (ielep,file=elename,status='old')
 701       call getenv('SIDEPAR',sidename)
 702       open (isidep,file=sidename,status='old')
 703       call getenv('SIDEP',sidepname)
 704       open (isidep1,file=sidepname,status="old")
 705       call getenv('SCCORPAR',sccorname)
 706       open (isccor,file=sccorname,status="old")
 707 #ifndef OLDSCP
 708 C
 709 C 8/9/01 In the newest version SCp interaction constants are read from a file
 710 C Use -DOLDSCP to use hard-coded constants instead.
 711 C
 712       call getenv('SCPPAR',scpname)
 713       open (iscpp,file=scpname,status='old')
 714 #endif
 715       return
 716       end
 717       subroutine read_dist_constr
 718       implicit real*8 (a-h,o-z)
 719       include 'DIMENSIONS'
 720 #ifdef MPI
 721       include 'mpif.h'
 722 #endif
 723       include 'COMMON.CONTROL'
 724       include 'COMMON.CHAIN'
 725       include 'COMMON.IOUNITS'
 726       include 'COMMON.SBRIDGE'
 727       integer ifrag_(2,100),ipair_(2,100)
 728       double precision wfrag_(100),wpair_(100)
 729       character*500 controlcard
 730       logical lprn /.true./
 731       write (iout,*) "Calling read_dist_constr"
 732 C      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 733 C      call flush(iout)
 734       write(iout,*) "TU sie wywalam?"
 735       call card_concat(controlcard)
 736       write (iout,*) controlcard
 737       call flush(iout)
 738       call readi(controlcard,"NFRAG",nfrag_,0)
 739       call readi(controlcard,"NPAIR",npair_,0)
 740       call readi(controlcard,"NDIST",ndist_,0)
 741       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 742       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 743       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 744       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 745       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 746       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 747       write (iout,*) "IFRAG"
 748       do i=1,nfrag_
 749         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 750       enddo
 751       write (iout,*) "IPAIR"
 752       do i=1,npair_
 753         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 754       enddo
 755       call flush(iout)
 756       do i=1,nfrag_
 757         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 758         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 759      &    ifrag_(2,i)=nstart_sup+nsup-1
 760 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 761         call flush(iout)
 762         if (wfrag_(i).gt.0.0d0) then
 763         do j=ifrag_(1,i),ifrag_(2,i)-1
 764           do k=j+1,ifrag_(2,i)
 765             write (iout,*) "j",j," k",k
 766             ddjk=dist(j,k)
 767             if (constr_dist.eq.1) then
 768               nhpb=nhpb+1
 769               ihpb(nhpb)=j
 770               jhpb(nhpb)=k
 771               dhpb(nhpb)=ddjk
 772               forcon(nhpb)=wfrag_(i)
 773             else if (constr_dist.eq.2) then
 774               if (ddjk.le.dist_cut) then
 775                 nhpb=nhpb+1
 776                 ihpb(nhpb)=j
 777                 jhpb(nhpb)=k
 778                 dhpb(nhpb)=ddjk
 779                 forcon(nhpb)=wfrag_(i)
 780               endif
 781             else
 782               nhpb=nhpb+1
 783               ihpb(nhpb)=j
 784               jhpb(nhpb)=k
 785               dhpb(nhpb)=ddjk
 786               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 787             endif
 788             if (lprn)
 789      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 790      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 791           enddo
 792         enddo
 793         endif
 794       enddo
 795       do i=1,npair_
 796         if (wpair_(i).gt.0.0d0) then
 797         ii = ipair_(1,i)
 798         jj = ipair_(2,i)
 799         if (ii.gt.jj) then
 800           itemp=ii
 801           ii=jj
 802           jj=itemp
 803         endif
 804         do j=ifrag_(1,ii),ifrag_(2,ii)
 805           do k=ifrag_(1,jj),ifrag_(2,jj)
 806             nhpb=nhpb+1
 807             ihpb(nhpb)=j
 808             jhpb(nhpb)=k
 809             forcon(nhpb)=wpair_(i)
 810             dhpb(nhpb)=dist(j,k)
 811             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 812      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 813           enddo
 814         enddo
 815         endif
 816       enddo
 817       do i=1,ndist_
 818         if (constr_dist.eq.11) then
 819         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 820      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 821         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 822 C        write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 823 C     &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 824         else
 825         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 826         endif
 827         if (forcon(nhpb+1).gt.0.0d0) then
 828           nhpb=nhpb+1
 829           if (ibecarb(i).gt.0) then
 830             ihpb(i)=ihpb(i)+nres
 831             jhpb(i)=jhpb(i)+nres
 832           endif
 833           if (dhpb(nhpb).eq.0.0d0)
 834      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 835 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 836           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 837      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 838         endif
 839 C      endif
 840       enddo
 841       call hpb_partition
 842       call flush(iout)
 843       return
 844       end