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