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