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