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