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
  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
 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 C 12/1/95 Added weight for the multi-body term WCORR
 130       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 131       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 132       weights(1)=wsc
 133       weights(2)=wscp
 134       weights(3)=welec
 135       weights(4)=wcorr
 136       weights(5)=wcorr5
 137       weights(6)=wcorr6
 138       weights(7)=wel_loc
 139       weights(8)=wturn3
 140       weights(9)=wturn4
 141       weights(10)=wturn6
 142       weights(11)=wang
 143       weights(12)=wscloc
 144       weights(13)=wtor
 145       weights(14)=wtor_d
 146       weights(15)=wstrain
 147       weights(16)=wvdwpp
 148       weights(17)=wbond
 149       weights(18)=scal14
 150       weights(19)=wsccor
 151       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 152      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 153      &  wturn4,wturn6,wsccor
 154    10 format (/'Energy-term weights (unscaled):'//
 155      & 'WSCC=   ',f10.6,' (SC-SC)'/
 156      & 'WSCP=   ',f10.6,' (SC-p)'/
 157      & 'WELEC=  ',f10.6,' (p-p electr)'/
 158      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 159      & 'WBOND=  ',f10.6,' (stretching)'/
 160      & 'WANG=   ',f10.6,' (bending)'/
 161      & 'WSCLOC= ',f10.6,' (SC local)'/
 162      & 'WTOR=   ',f10.6,' (torsional)'/
 163      & 'WTORD=  ',f10.6,' (double torsional)'/
 164      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 165      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 166      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 167      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 168      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 169      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 170      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 171      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 172      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 173
 174       if (wcorr4.gt.0.0d0) then
 175         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 176      &   'between contact pairs of peptide groups'
 177         write (iout,'(2(a,f5.3/))')
 178      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 179      &  'Range of quenching the correlation terms:',2*delt_corr
 180       else if (wcorr.gt.0.0d0) then
 181         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 182      &   'between contact pairs of peptide groups'
 183       endif
 184       write (iout,'(a,f8.3)')
 185      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 186       write (iout,'(a,f8.3)')
 187      &  'General scaling factor of SC-p interactions:',scalscp
 188       r0_corr=cutoff_corr-delt_corr
 189       do i=1,20
 190         aad(i,1)=scalscp*aad(i,1)
 191         aad(i,2)=scalscp*aad(i,2)
 192         bad(i,1)=scalscp*bad(i,1)
 193         bad(i,2)=scalscp*bad(i,2)
 194       enddo
 195
 196       call flush(iout)
 197       print *,'indpdb=',indpdb,' pdbref=',pdbref
 198
 199 C Read sequence if not taken from the pdb file.
 200       if (iscode.gt.0) then
 201         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 202       else
 203         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 204       endif
 205 C Convert sequence to numeric code
 206       do i=1,nres
 207         itype(i)=rescode(i,sequence(i),iscode)
 208       enddo
 209       print *,nres
 210       print '(20i4)',(itype(i),i=1,nres)
 211
 212       do i=1,nres
 213 #ifdef PROCOR
 214         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 215 #else
 216         if (itype(i).eq.21) then
 217 #endif
 218           itel(i)=0
 219 #ifdef PROCOR
 220         else if (itype(i+1).ne.20) then
 221 #else
 222         else if (itype(i).ne.20) then
 223 #endif
 224           itel(i)=1
 225         else
 226           itel(i)=2
 227         endif
 228       enddo
 229       write (iout,*) "ITEL"
 230       do i=1,nres-1
 231         write (iout,*) i,itype(i),itel(i)
 232       enddo
 233
 234       print *,'Call Read_Bridge.'
 235       call read_bridge
 236       nnt=1
 237       nct=nres
 238       print *,'NNT=',NNT,' NCT=',NCT
 239       if (itype(1).eq.21) nnt=2
 240       if (itype(nres).eq.21) nct=nct-1
 241       if (nstart.lt.nnt) nstart=nnt
 242       if (nend.gt.nct .or. nend.eq.0) nend=nct
 243       write (iout,*) "nstart",nstart," nend",nend
 244       nres0=nres
 245 c      if (pdbref) then
 246 c        read(inp,'(a)') pdbfile
 247 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 248 c        open(ipdbin,file=pdbfile,status='old',err=33)
 249 c        goto 34
 250 c  33    write (iout,'(a)') 'Error opening PDB file.'
 251 c        stop
 252 c  34    continue
 253 c        print *,'Begin reading pdb data'
 254 c        call readpdb
 255 c        print *,'Finished reading pdb data'
 256 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 257 c        do i=1,nres
 258 c          itype_pdb(i)=itype(i)
 259 c        enddo
 260 c        close (ipdbin)
 261 c        write (iout,'(a,i3)') 'nsup=',nsup
 262 c        nstart_seq=nnt
 263 c        if (nsup.le.(nct-nnt+1)) then
 264 c          do i=0,nct-nnt+1-nsup
 265 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 266 c              nstart_seq=nnt+i
 267 c              goto 111
 268 c            endif
 269 c          enddo
 270 c          write (iout,'(a)')
 271 c     &            'Error - sequences to be superposed do not match.'
 272 c          stop
 273 c        else
 274 c          do i=0,nsup-(nct-nnt+1)
 275 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 276 c     &      then
 277 c              nstart_sup=nstart_sup+i
 278 c              nsup=nct-nnt+1
 279 c              goto 111
 280 c            endif
 281 c          enddo
 282 c          write (iout,'(a)')
 283 c     &            'Error - sequences to be superposed do not match.'
 284 c        endif
 285 c  111   continue
 286 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 287 c     &                 ' nstart_seq=',nstart_seq
 288 c      endif
 289       call init_int_table
 290       call setup_var
 291       write (iout,*) "molread: REFSTR",refstr
 292       if (refstr) then
 293         if (.not.pdbref) then
 294           call read_angles(inp,*38)
 295           goto 39
 296    38     write (iout,'(a)') 'Error reading reference structure.'
 297 #ifdef MPL
 298           call mp_stopall(Error_Msg)
 299 #else
 300           stop 'Error reading reference structure'
 301 #endif
 302    39     call chainbuild
 303           nstart_sup=nnt
 304           nstart_seq=nnt
 305           nsup=nct-nnt+1
 306           do i=1,2*nres
 307             do j=1,3
 308               cref(j,i)=c(j,i)
 309             enddo
 310           enddo
 311         endif
 312         call contact(.true.,ncont_ref,icont_ref)
 313       endif
 314       return
 315       end
 316 c-----------------------------------------------------------------------------
 317       logical function seq_comp(itypea,itypeb,length)
 318       implicit none
 319       integer length,itypea(length),itypeb(length)
 320       integer i
 321       do i=1,length
 322         if (itypea(i).ne.itypeb(i)) then
 323           seq_comp=.false.
 324           return
 325         endif
 326       enddo
 327       seq_comp=.true.
 328       return
 329       end
 330 c-----------------------------------------------------------------------------
 331       subroutine read_bridge
 332 C Read information about disulfide bridges.
 333       implicit none
 334       include 'DIMENSIONS'
 335       include 'COMMON.IOUNITS'
 336       include 'COMMON.GEO'
 337       include 'COMMON.VAR'
 338       include 'COMMON.INTERACT'
 339       include 'COMMON.LOCAL'
 340       include 'COMMON.NAMES'
 341       include 'COMMON.CHAIN'
 342       include 'COMMON.FFIELD'
 343       include 'COMMON.SBRIDGE'
 344       include 'COMMON.HEADER'
 345       include 'COMMON.CONTROL'
 346       include 'COMMON.TIME1'
 347 #ifdef MPL
 348       include 'COMMON.INFO'
 349 #endif
 350       integer i,j
 351 C Read bridging residues.
 352       read (inp,*) ns,(iss(i),i=1,ns)
 353       print *,'ns=',ns
 354 C Check whether the specified bridging residues are cystines.
 355       do i=1,ns
 356         if (itype(iss(i)).ne.1) then
 357           write (iout,'(2a,i3,a)')
 358      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 359      &   ' can form a disulfide bridge?!!!'
 360           write (*,'(2a,i3,a)')
 361      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 362      &   ' can form a disulfide bridge?!!!'
 363 #ifdef MPL
 364          call mp_stopall(error_msg)
 365 #else
 366          stop
 367 #endif
 368         endif
 369       enddo
 370 C Read preformed bridges.
 371       if (ns.gt.0) then
 372       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 373       if (nss.gt.0) then
 374         nhpb=nss
 375 C Check if the residues involved in bridges are in the specified list of
 376 C bridging residues.
 377         do i=1,nss
 378           do j=1,i-1
 379             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 380      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 381               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 382      &      ' contains residues present in other pairs.'
 383               write (*,'(a,i3,a)') 'Disulfide pair',i,
 384      &      ' contains residues present in other pairs.'
 385 #ifdef MPL
 386               call mp_stopall(error_msg)
 387 #else
 388               stop
 389 #endif
 390             endif
 391           enddo
 392           do j=1,ns
 393             if (ihpb(i).eq.iss(j)) goto 10
 394           enddo
 395           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 396    10     continue
 397           do j=1,ns
 398             if (jhpb(i).eq.iss(j)) goto 20
 399           enddo
 400           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 401    20     continue
 402           dhpb(i)=dbr
 403           forcon(i)=fbr
 404         enddo
 405         do i=1,nss
 406           ihpb(i)=ihpb(i)+nres
 407           jhpb(i)=jhpb(i)+nres
 408         enddo
 409       endif
 410       endif
 411       return
 412       end
 413 c----------------------------------------------------------------------------
 414       subroutine read_angles(kanal,*)
 415       implicit none
 416       include 'DIMENSIONS'
 417       include 'COMMON.GEO'
 418       include 'COMMON.VAR'
 419       include 'COMMON.CHAIN'
 420       include 'COMMON.IOUNITS'
 421       integer i,kanal
 422       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 423       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 424       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 425       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 426       do i=1,nres
 427         theta(i)=deg2rad*theta(i)
 428         phi(i)=deg2rad*phi(i)
 429         alph(i)=deg2rad*alph(i)
 430         omeg(i)=deg2rad*omeg(i)
 431       enddo
 432       return
 433    10 return1
 434       end
 435 c----------------------------------------------------------------------------
 436       subroutine reada(rekord,lancuch,wartosc,default)
 437       implicit none
 438       character*(*) rekord,lancuch
 439       double precision wartosc,default
 440       integer ilen,iread
 441       external ilen
 442       iread=index(rekord,lancuch)
 443       if (iread.eq.0) then
 444         wartosc=default
 445         return
 446       endif
 447       iread=iread+ilen(lancuch)+1
 448       read (rekord(iread:),*) wartosc
 449       return
 450       end
 451 c----------------------------------------------------------------------------
 452       subroutine multreada(rekord,lancuch,tablica,dim,default)
 453       implicit none
 454       integer dim,i
 455       double precision tablica(dim),default
 456       character*(*) rekord,lancuch
 457       integer ilen,iread
 458       external ilen
 459       do i=1,dim
 460         tablica(i)=default
 461       enddo
 462       iread=index(rekord,lancuch)
 463       if (iread.eq.0) return
 464       iread=iread+ilen(lancuch)+1
 465       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 466    10 return
 467       end
 468 c----------------------------------------------------------------------------
 469       subroutine readi(rekord,lancuch,wartosc,default)
 470       implicit none
 471       character*(*) rekord,lancuch
 472       integer wartosc,default
 473       integer ilen,iread
 474       external ilen
 475       iread=index(rekord,lancuch)
 476       if (iread.eq.0) then
 477         wartosc=default
 478         return
 479       endif
 480       iread=iread+ilen(lancuch)+1
 481       read (rekord(iread:),*) wartosc
 482       return
 483       end
 484 c----------------------------------------------------------------------------
 485       subroutine card_concat(card)
 486       include 'DIMENSIONS'
 487       include 'COMMON.IOUNITS'
 488       character*(*) card
 489       character*80 karta,ucase
 490       external ilen
 491       read (inp,'(a)') karta
 492       karta=ucase(karta)
 493       card=' '
 494       do while (karta(80:80).eq.'&')
 495         card=card(:ilen(card)+1)//karta(:79)
 496         read (inp,'(a)') karta
 497         karta=ucase(karta)
 498       enddo
 499       card=card(:ilen(card)+1)//karta
 500       return
 501       end
 502 c----------------------------------------------------------------------------
 503       subroutine openunits
 504       implicit none
 505       include 'DIMENSIONS'
 506 #ifdef MPI
 507       include "mpif.h"
 508       character*3 liczba
 509       include "COMMON.MPI"
 510 #endif
 511       include 'COMMON.IOUNITS'
 512       include 'COMMON.CONTROL'
 513       integer lenpre,lenpot,ilen
 514       external ilen
 515       character*16 cformat,cprint
 516       character*16 ucase
 517       integer lenint,lenout
 518       call getenv('INPUT',prefix)
 519       call getenv('OUTPUT',prefout)
 520       call getenv('INTIN',prefintin)
 521       call getenv('COORD',cformat)
 522       call getenv('PRINTCOOR',cprint)
 523       call getenv('SCRATCHDIR',scratchdir)
 524       from_bx=.true.
 525       from_cx=.false.
 526       if (index(ucase(cformat),'CX').gt.0) then
 527         from_cx=.true.
 528         from_bx=.false.
 529       endif
 530       from_cart=.true.
 531       lenpre=ilen(prefix)
 532       lenout=ilen(prefout)
 533       lenint=ilen(prefintin)
 534 C Get the names and open the input files
 535       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 536 #ifdef MPI
 537       write (liczba,'(bz,i3.3)') me
 538       outname=prefout(:lenout)//'_clust.out_'//liczba
 539 #else
 540       outname=prefout(:lenout)//'_clust.out'
 541 #endif
 542       if (from_bx) then
 543         intinname=prefintin(:lenint)//'.bx'
 544       else if (from_cx) then
 545         intinname=prefintin(:lenint)//'.cx'
 546       else
 547         intinname=prefintin(:lenint)//'.int'
 548       endif
 549       rmsname=prefintin(:lenint)//'.rms'
 550       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 551      &       status='unknown')
 552       open (jrms,file=rmsname,status='unknown')
 553       open(iout,file=outname,status='unknown')
 554 C Get parameter filenames and open the parameter files.
 555       call getenv('BONDPAR',bondname)
 556       open (ibond,file=bondname,status='old')
 557       call getenv('THETPAR',thetname)
 558       open (ithep,file=thetname,status='old')
 559       call getenv('ROTPAR',rotname)
 560       open (irotam,file=rotname,status='old')
 561       call getenv('TORPAR',torname)
 562       open (itorp,file=torname,status='old')
 563       call getenv('TORDPAR',tordname)
 564       open (itordp,file=tordname,status='old')
 565       call getenv('FOURIER',fouriername)
 566       open (ifourier,file=fouriername,status='old')
 567       call getenv('ELEPAR',elename)
 568       open (ielep,file=elename,status='old')
 569       call getenv('SIDEPAR',sidename)
 570       open (isidep,file=sidename,status='old')
 571       call getenv('SIDEP',sidepname)
 572       open (isidep1,file=sidepname,status="old")
 573       call getenv('SCCORPAR',sccorname)
 574       open (isccor,file=sccorname,status="old")
 575 #ifndef OLDSCP
 576 C
 577 C 8/9/01 In the newest version SCp interaction constants are read from a file
 578 C Use -DOLDSCP to use hard-coded constants instead.
 579 C
 580       call getenv('SCPPAR',scpname)
 581       open (iscpp,file=scpname,status='old')
 582 #endif
 583       return
 584       end