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