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