source/unres/src_MD/parmread.F

   1       subroutine parmread
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6 C Important! Energy-term weights ARE NOT read here; they are read from the
   7 C main input file instead, because NO defaults have yet been set for these
   8 C parameters.
   9 C
  10       implicit real*8 (a-h,o-z)
  11       include 'DIMENSIONS'
  12 #ifdef MPI
  13       include "mpif.h"
  14       integer IERROR
  15 #endif
  16       include 'COMMON.IOUNITS'
  17       include 'COMMON.CHAIN'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.GEO'
  20       include 'COMMON.LOCAL'
  21       include 'COMMON.TORSION'
  22       include 'COMMON.SCCOR'
  23       include 'COMMON.SCROT'
  24       include 'COMMON.FFIELD'
  25       include 'COMMON.NAMES'
  26       include 'COMMON.SBRIDGE'
  27       include 'COMMON.MD'
  28       include 'COMMON.SETUP'
  29       character*1 t1,t2,t3
  30       character*1 onelett(4) /"G","A","P","D"/
  31       logical lprint,LaTeX
  32       dimension blower(3,3,maxlob)
  33       dimension b(13)
  34       character*3 lancuch,ucase
  35 C
  36 C For printing parameters after they are read set the following in the UNRES
  37 C C-shell script:
  38 C
  39 C setenv PRINT_PARM YES
  40 C
  41 C To print parameters in LaTeX format rather than as ASCII tables:
  42 C
  43 C setenv LATEX YES
  44 C
  45       call getenv_loc("PRINT_PARM",lancuch)
  46       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  47       call getenv_loc("LATEX",lancuch)
  48       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  49 C
  50       dwa16=2.0d0**(1.0d0/6.0d0)
  51       itypro=20
  52 C Assign virtual-bond length
  53       vbl=3.8D0
  54       vblinv=1.0D0/vbl
  55       vblinv2=vblinv*vblinv
  56 c
  57 c Read the virtual-bond parameters, masses, and moments of inertia
  58 c and Stokes' radii of the peptide group and side chains
  59 c
  60 #ifdef CRYST_BOND
  61       read (ibond,*) vbldp0,akp,mp,ip,pstok
  62       do i=1,ntyp
  63         nbondterm(i)=1
  64         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
  65         dsc(i) = vbldsc0(1,i)
  66         if (i.eq.10) then
  67           dsc_inv(i)=0.0D0
  68         else
  69           dsc_inv(i)=1.0D0/dsc(i)
  70         endif
  71       enddo
  72 #else
  73       read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
  74       do i=1,ntyp
  75         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
  76      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
  77         dsc(i) = vbldsc0(1,i)
  78         if (i.eq.10) then
  79           dsc_inv(i)=0.0D0
  80         else
  81           dsc_inv(i)=1.0D0/dsc(i)
  82         endif
  83       enddo
  84 #endif
  85       if (lprint) then
  86         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
  87         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
  88      &   'inertia','Pstok'
  89         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
  90         do i=1,ntyp
  91           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
  92      &      vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
  93           do j=2,nbondterm(i)
  94             write (iout,'(13x,3f10.5)')
  95      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
  96           enddo
  97         enddo
  98       endif
  99 #ifdef CRYST_THETA
 100 C
 101 C Read the parameters of the probability distribution/energy expression
 102 C of the virtual-bond valence angles theta
 103 C
 104       do i=1,ntyp
 105         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
 106      &    (bthet(j,i),j=1,2)
 107         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 108         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 109         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 110         sigc0(i)=sigc0(i)**2
 111       enddo
 112       close (ithep)
 113       if (lprint) then
 114       if (.not.LaTeX) then
 115         write (iout,'(a)')
 116      &    'Parameters of the virtual-bond valence angles:'
 117         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 118      & '    ATHETA0   ','         A1   ','        A2    ',
 119      & '        B1    ','         B2   '
 120         do i=1,ntyp
 121           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 122      &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 123         enddo
 124         write (iout,'(/a/9x,5a/79(1h-))')
 125      & 'Parameters of the expression for sigma(theta_c):',
 126      & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 127      & '     ALPH3    ','    SIGMA0C   '
 128         do i=1,ntyp
 129           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 130      &      (polthet(j,i),j=0,3),sigc0(i)
 131         enddo
 132         write (iout,'(/a/9x,5a/79(1h-))')
 133      & 'Parameters of the second gaussian:',
 134      & '    THETA0    ','     SIGMA0   ','        G1    ',
 135      & '        G2    ','         G3   '
 136         do i=1,ntyp
 137           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 138      &       sig0(i),(gthet(j,i),j=1,3)
 139         enddo
 140        else
 141         write (iout,'(a)')
 142      &    'Parameters of the virtual-bond valence angles:'
 143         write (iout,'(/a/9x,5a/79(1h-))')
 144      & 'Coefficients of expansion',
 145      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 146      & '   b1*10^1    ','    b2*10^1   '
 147         do i=1,ntyp
 148           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 149      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
 150         enddo
 151         write (iout,'(/a/9x,5a/79(1h-))')
 152      & 'Parameters of the expression for sigma(theta_c):',
 153      & ' alpha0       ','  alph1       ',' alph2        ',
 154      & ' alhp3        ','   sigma0c    '
 155         do i=1,ntyp
 156           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 157      &      (polthet(j,i),j=0,3),sigc0(i)
 158         enddo
 159         write (iout,'(/a/9x,5a/79(1h-))')
 160      & 'Parameters of the second gaussian:',
 161      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 162      & '        G2    ','   G3*10^1    '
 163         do i=1,ntyp
 164           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 165      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 166         enddo
 167       endif
 168       endif
 169 #else
 170 C
 171 C Read the parameters of Utheta determined from ab initio surfaces
 172 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 173 C
 174       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
 175      &  ntheterm3,nsingle,ndouble
 176       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 177       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
 178       do i=1,maxthetyp
 179         do j=1,maxthetyp
 180           do k=1,maxthetyp
 181             aa0thet(i,j,k)=0.0d0
 182             do l=1,ntheterm
 183               aathet(l,i,j,k)=0.0d0
 184             enddo
 185             do l=1,ntheterm2
 186               do m=1,nsingle
 187                 bbthet(m,l,i,j,k)=0.0d0
 188                 ccthet(m,l,i,j,k)=0.0d0
 189                 ddthet(m,l,i,j,k)=0.0d0
 190                 eethet(m,l,i,j,k)=0.0d0
 191               enddo
 192             enddo
 193             do l=1,ntheterm3
 194               do m=1,ndouble
 195                 do mm=1,ndouble
 196                  ffthet(mm,m,l,i,j,k)=0.0d0
 197                  ggthet(mm,m,l,i,j,k)=0.0d0
 198                 enddo
 199               enddo
 200             enddo
 201           enddo
 202         enddo
 203       enddo
 204       do i=1,nthetyp
 205         do j=1,nthetyp
 206           do k=1,nthetyp
 207             read (ithep,'(3a)',end=111,err=111) res1,res2,res3
 208             read (ithep,*,end=111,err=111) aa0thet(i,j,k)
 209             read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
 210             read (ithep,*,end=111,err=111)
 211      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
 212      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
 213      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
 214      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
 215             read (ithep,*,end=111,err=111)
 216      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
 217      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
 218      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 219           enddo
 220         enddo
 221       enddo
 222 C
 223 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 224 C coefficients of theta-and-gamma-dependent terms are zero.
 225 C
 226       do i=1,nthetyp
 227         do j=1,nthetyp
 228           do l=1,ntheterm
 229             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
 230             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
 231           enddo
 232           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
 233           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
 234         enddo
 235         do l=1,ntheterm
 236           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
 237         enddo
 238         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
 239       enddo
 240 C
 241 C Control printout of the coefficients of virtual-bond-angle potentials
 242 C
 243       if (lprint) then
 244         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 245         do i=1,nthetyp+1
 246           do j=1,nthetyp+1
 247             do k=1,nthetyp+1
 248               write (iout,'(//4a)')
 249      &         'Type ',onelett(i),onelett(j),onelett(k)
 250               write (iout,'(//a,10x,a)') " l","a[l]"
 251               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
 252               write (iout,'(i2,1pe15.5)')
 253      &           (l,aathet(l,i,j,k),l=1,ntheterm)
 254             do l=1,ntheterm2
 255               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
 256      &          "b",l,"c",l,"d",l,"e",l
 257               do m=1,nsingle
 258                 write (iout,'(i2,4(1pe15.5))') m,
 259      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
 260      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
 261               enddo
 262             enddo
 263             do l=1,ntheterm3
 264               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
 265      &          "f+",l,"f-",l,"g+",l,"g-",l
 266               do m=2,ndouble
 267                 do n=1,m-1
 268                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 269      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
 270      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
 271                 enddo
 272               enddo
 273             enddo
 274           enddo
 275         enddo
 276       enddo
 277       call flush(iout)
 278       endif
 279       write (2,*) "Start reading THETA_PDB"
 280       do i=1,ntyp
 281         read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
 282      &    (bthet(j,i),j=1,2)
 283         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
 284         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
 285         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 286         sigc0(i)=sigc0(i)**2
 287       enddo
 288       write (2,*) "End reading THETA_PDB"
 289       close (ithep_pdb)
 290 #endif
 291       close(ithep)
 292 #ifdef CRYST_SC
 293 C
 294 C Read the parameters of the probability distribution/energy expression
 295 C of the side chains.
 296 C
 297       do i=1,ntyp
 298         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 299         if (i.eq.10) then
 300           dsc_inv(i)=0.0D0
 301         else
 302           dsc_inv(i)=1.0D0/dsc(i)
 303         endif
 304         if (i.ne.10) then
 305         do j=1,nlob(i)
 306           do k=1,3
 307             do l=1,3
 308               blower(l,k,j)=0.0D0
 309             enddo
 310           enddo
 311         enddo
 312         bsc(1,i)=0.0D0
 313         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 314      &    ((blower(k,l,1),l=1,k),k=1,3)
 315         censc(1,1,-i)=censc(1,1,i)
 316         censc(2,1,-i)=censc(2,1,i)
 317         censc(3,1,-i)=-censc(3,1,i)
 318
 319         do j=2,nlob(i)
 320           read (irotam,*,end=112,err=112) bsc(j,i)
 321           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 322      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 323         censc(1,j,-i)=censc(1,j,i)
 324         censc(2,j,-i)=censc(2,j,i)
 325         censc(3,j,-i)=-censc(3,j,i)
 326 C BSC is amplitude of Gaussian
 327         enddo
 328         do j=1,nlob(i)
 329           do k=1,3
 330             do l=1,k
 331               akl=0.0D0
 332               do m=1,3
 333                 akl=akl+blower(k,m,j)*blower(l,m,j)
 334               enddo
 335               gaussc(k,l,j,i)=akl
 336               gaussc(l,k,j,i)=akl
 337               if (((k.eq.3).and.(l.ne.3))
 338      &        .or.((l.eq.3).and.(k.ne.3))) then
 339                 gaussc(k,l,j,-i)=-akl
 340                 gaussc(l,k,j,-i)=-akl
 341               else
 342                 gaussc(k,l,j,-i)=akl
 343                 gaussc(l,k,j,-i)=akl
 344               endif
 345             enddo
 346           enddo
 347         enddo
 348         endif
 349       enddo
 350       close (irotam)
 351       if (lprint) then
 352         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 353         do i=1,ntyp
 354           nlobi=nlob(i)
 355           if (nlobi.gt.0) then
 356             if (LaTeX) then
 357               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 358      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 359                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 360      &                             'log h',(bsc(j,i),j=1,nlobi)
 361                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 362      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 363               do k=1,3
 364                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 365      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 366               enddo
 367             else
 368               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 369               write (iout,'(a,f10.4,4(16x,f10.4))')
 370      &                             'Center  ',(bsc(j,i),j=1,nlobi)
 371               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
 372      &           j=1,nlobi)
 373               write (iout,'(a)')
 374             endif
 375           endif
 376         enddo
 377       endif
 378 #else
 379 C
 380 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 381 C added by Urszula Kozlowska 07/11/2007
 382 C
 383       do i=1,ntyp
 384         read (irotam,*,end=112,err=112)
 385        if (i.eq.10) then
 386          read (irotam,*,end=112,err=112)
 387        else
 388          do j=1,65
 389            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 390          enddo
 391        endif
 392       enddo
 393 C
 394 C Read the parameters of the probability distribution/energy expression
 395 C of the side chains.
 396 C
 397       do i=1,ntyp
 398         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 399         if (i.eq.10) then
 400           dsc_inv(i)=0.0D0
 401         else
 402           dsc_inv(i)=1.0D0/dsc(i)
 403         endif
 404         if (i.ne.10) then
 405         do j=1,nlob(i)
 406           do k=1,3
 407             do l=1,3
 408               blower(l,k,j)=0.0D0
 409             enddo
 410           enddo
 411         enddo
 412         bsc(1,i)=0.0D0
 413         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
 414      &    ((blower(k,l,1),l=1,k),k=1,3)
 415         do j=2,nlob(i)
 416           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
 417           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
 418      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 419         enddo
 420         do j=1,nlob(i)
 421           do k=1,3
 422             do l=1,k
 423               akl=0.0D0
 424               do m=1,3
 425                 akl=akl+blower(k,m,j)*blower(l,m,j)
 426               enddo
 427               gaussc(k,l,j,i)=akl
 428               gaussc(l,k,j,i)=akl
 429             enddo
 430           enddo
 431         enddo
 432         endif
 433       enddo
 434       close (irotam_pdb)
 435 #endif
 436       close(irotam)
 437
 438 #ifdef CRYST_TOR
 439 C
 440 C Read torsional parameters in old format
 441 C
 442       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 443       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 444       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 445       do i=1,ntortyp
 446         do j=1,ntortyp
 447           read (itorp,'(a)')
 448           do k=1,nterm_old
 449             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 450           enddo
 451         enddo
 452       enddo
 453       close (itorp)
 454       if (lprint) then
 455         write (iout,'(/a/)') 'Torsional constants:'
 456         do i=1,ntortyp
 457           do j=1,ntortyp
 458             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 459             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 460           enddo
 461         enddo
 462       endif
 463 #else
 464 C
 465 C Read torsional parameters
 466 C
 467       read (itorp,*,end=113,err=113) ntortyp
 468       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 469       do iblock=1,2
 470       do i=-ntyp,-1
 471        itortyp(i)=-itortyp(-i)
 472       enddo
 473 c      write (iout,*) 'ntortyp',ntortyp
 474       do i=0,ntortyp-1
 475         do j=-ntortyp,ntortyp
 476           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 477      &          nlor(i,j,iblock)
 478           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 479           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 480           v0ij=0.0d0
 481           si=-1.0d0
 482           do k=1,nterm(i,j,iblock)
 483             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 484      &      v2(k,i,j,iblock)
 485             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 486             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 487             v0ij=v0ij+si*v1(k,i,j)
 488             si=-si
 489           enddo
 490           do k=1,nlor(i,j,iblock)
 491             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 492      &        vlor2(k,i,j),vlor3(k,i,j)
 493             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 494           enddo
 495           v0(i,j,iblock)=v0ij
 496           v0(-i,-j,iblock)=v0ij
 497         enddo
 498       enddo
 499       enddo
 500       close (itorp)
 501       if (lprint) then
 502         write (iout,'(/a/)') 'Torsional constants:'
 503         do i=1,ntortyp
 504           do j=1,ntortyp
 505             write (iout,*) 'ityp',i,' jtyp',j
 506             write (iout,*) 'Fourier constants'
 507             do k=1,nterm(i,j)
 508               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 509             enddo
 510             write (iout,*) 'Lorenz constants'
 511             do k=1,nlor(i,j)
 512               write (iout,'(3(1pe15.5))')
 513      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 514             enddo
 515           enddo
 516         enddo
 517       endif
 518 C
 519 C 6/23/01 Read parameters for double torsionals
 520 C
 521       do i=1,ntortyp
 522         do j=1,ntortyp
 523           do k=1,ntortyp
 524             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 525             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 526      &        .or. t3.ne.onelett(k)) then
 527               write (iout,*) "Error in double torsional parameter file",
 528      &         i,j,k,t1,t2,t3
 529 #ifdef MPI
 530               call MPI_Finalize(Ierror)
 531 #endif
 532                stop "Error in double torsional parameter file"
 533             endif
 534             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
 535      &         ntermd_2(i,j,k)
 536             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
 537      &         ntermd_1(i,j,k))
 538             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
 539      &         ntermd_1(i,j,k))
 540             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
 541      &         ntermd_1(i,j,k))
 542             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
 543      &         ntermd_1(i,j,k))
 544             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
 545      &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
 546      &         m=1,l-1),l=1,ntermd_2(i,j,k))
 547           enddo
 548         enddo
 549       enddo
 550       if (lprint) then
 551       write (iout,*)
 552       write (iout,*) 'Constants for double torsionals'
 553       do i=1,ntortyp
 554         do j=1,ntortyp
 555           do k=1,ntortyp
 556             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 557      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 558             write (iout,*)
 559             write (iout,*) 'Single angles:'
 560             do l=1,ntermd_1(i,j,k)
 561               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 562      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 563      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 564             enddo
 565             write (iout,*)
 566             write (iout,*) 'Pairs of angles:'
 567             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 568             do l=1,ntermd_2(i,j,k)
 569               write (iout,'(i5,20f10.5)')
 570      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 571             enddo
 572             write (iout,*)
 573             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 574             do l=1,ntermd_2(i,j,k)
 575               write (iout,'(i5,20f10.5)')
 576      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 577             enddo
 578             write (iout,*)
 579           enddo
 580         enddo
 581       enddo
 582       endif
 583 #endif
 584 C Read of Side-chain backbone correlation parameters
 585 C Modified 11 May 2012 by Adasko
 586 CCC
 587 C
 588       read (isccor,*,end=113,err=113) nsccortyp
 589 #ifdef SCCORPDB
 590       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
 591 #else
 592       read (isccor,*,end=113,err=113) (isccortyp(i),i=-ntyp,ntyp)
 593 #endif
 594 c      write (iout,*) 'ntortyp',ntortyp
 595       maxinter=3
 596 cc maxinter is maximum interaction sites
 597       do l=1,maxinter
 598       do i=1,nsccortyp
 599         do j=1,nsccortyp
 600           read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
 601           v0ijsccor=0.0d0
 602           si=-1.0d0
 603
 604           do k=1,nterm_sccor(i,j)
 605             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
 606      &    ,v2sccor(k,l,i,j)
 607             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 608             si=-si
 609           enddo
 610           do k=1,nlor_sccor(i,j)
 611             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 612      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 613             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 614      &(1+vlor3sccor(k,i,j)**2)
 615           enddo
 616           v0sccor(i,j)=v0ijsccor
 617         enddo
 618       enddo
 619       enddo
 620       close (isccor)
 621
 622       if (lprint) then
 623         write (iout,'(/a/)') 'Torsional constants:'
 624         do i=1,nsccortyp
 625           do j=1,nsccortyp
 626             write (iout,*) 'ityp',i,' jtyp',j
 627             write (iout,*) 'Fourier constants'
 628             do k=1,nterm_sccor(i,j)
 629       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 630             enddo
 631             write (iout,*) 'Lorenz constants'
 632             do k=1,nlor_sccor(i,j)
 633               write (iout,'(3(1pe15.5))')
 634      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 635             enddo
 636           enddo
 637         enddo
 638       endif
 639 C
 640 C
 641 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 642 C         interaction energy of the Gly, Ala, and Pro prototypes.
 643 C
 644       if (lprint) then
 645         write (iout,*)
 646         write (iout,*) "Coefficients of the cumulants"
 647       endif
 648       read (ifourier,*) nloctyp
 649       do i=0,nloctyp-1
 650         read (ifourier,*,end=115,err=115)
 651         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 652         if (lprint) then
 653         write (iout,*) 'Type',i
 654         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 655         endif
 656         B1(1,i)  = b(3)
 657         B1(2,i)  = b(5)
 658         B1(1,-i) = b(3)
 659         B1(2,-i) = -b(5)
 660 c        b1(1,i)=0.0d0
 661 c        b1(2,i)=0.0d0
 662         B1tilde(1,i) = b(3)
 663         B1tilde(2,i) =-b(5)
 664         B1tilde(1,-i) =-b(3)
 665         B1tilde(2,-i) =b(5)
 666 c        b1tilde(1,i)=0.0d0
 667 c        b1tilde(2,i)=0.0d0
 668         B2(1,i)  = b(2)
 669         B2(2,i)  = b(4)
 670         B2(1,-i)  =b(2)
 671         B2(2,-i)  =-b(4)
 672
 673 c        b2(1,i)=0.0d0
 674 c        b2(2,i)=0.0d0
 675         CC(1,1,i)= b(7)
 676         CC(2,2,i)=-b(7)
 677         CC(2,1,i)= b(9)
 678         CC(1,2,i)= b(9)
 679         CC(1,1,-i)= b(7)
 680         CC(2,2,-i)=-b(7)
 681         CC(2,1,-i)=-b(9)
 682         CC(1,2,-i)=-b(9)
 683 c        CC(1,1,i)=0.0d0
 684 c        CC(2,2,i)=0.0d0
 685 c        CC(2,1,i)=0.0d0
 686 c        CC(1,2,i)=0.0d0
 687         Ctilde(1,1,i)=b(7)
 688         Ctilde(1,2,i)=b(9)
 689         Ctilde(2,1,i)=-b(9)
 690         Ctilde(2,2,i)=b(7)
 691         Ctilde(1,1,-i)=b(7)
 692         Ctilde(1,2,-i)=-b(9)
 693         Ctilde(2,1,-i)=b(9)
 694         Ctilde(2,2,-i)=b(7)
 695
 696 c        Ctilde(1,1,i)=0.0d0
 697 c        Ctilde(1,2,i)=0.0d0
 698 c        Ctilde(2,1,i)=0.0d0
 699 c        Ctilde(2,2,i)=0.0d0
 700         DD(1,1,i)= b(6)
 701         DD(2,2,i)=-b(6)
 702         DD(2,1,i)= b(8)
 703         DD(1,2,i)= b(8)
 704         DD(1,1,-i)= b(6)
 705         DD(2,2,-i)=-b(6)
 706         DD(2,1,-i)=-b(8)
 707         DD(1,2,-i)=-b(8)
 708 c        DD(1,1,i)=0.0d0
 709 c        DD(2,2,i)=0.0d0
 710 c        DD(2,1,i)=0.0d0
 711 c        DD(1,2,i)=0.0d0
 712         Dtilde(1,1,i)=b(6)
 713         Dtilde(1,2,i)=b(8)
 714         Dtilde(2,1,i)=-b(8)
 715         Dtilde(2,2,i)=b(6)
 716         Dtilde(1,1,-i)=b(6)
 717         Dtilde(1,2,-i)=-b(8)
 718         Dtilde(2,1,-i)=b(8)
 719         Dtilde(2,2,-i)=b(6)
 720
 721 c        Dtilde(1,1,i)=0.0d0
 722 c        Dtilde(1,2,i)=0.0d0
 723 c        Dtilde(2,1,i)=0.0d0
 724 c        Dtilde(2,2,i)=0.0d0
 725         EE(1,1,i)= b(10)+b(11)
 726         EE(2,2,i)=-b(10)+b(11)
 727         EE(2,1,i)= b(12)-b(13)
 728         EE(1,2,i)= b(12)+b(13)
 729         EE(1,1,-i)= b(10)+b(11)
 730         EE(2,2,-i)=-b(10)+b(11)
 731         EE(2,1,-i)=-b(12)+b(13)
 732         EE(1,2,-i)=-b(12)-b(13)
 733
 734 c        ee(1,1,i)=1.0d0
 735 c        ee(2,2,i)=1.0d0
 736 c        ee(2,1,i)=0.0d0
 737 c        ee(1,2,i)=0.0d0
 738 c        ee(2,1,i)=ee(1,2,i)
 739       enddo
 740       if (lprint) then
 741       do i=1,nloctyp
 742         write (iout,*) 'Type',i
 743         write (iout,*) 'B1'
 744         write(iout,*) B1(1,i),B1(2,i)
 745         write (iout,*) 'B2'
 746         write(iout,*) B2(1,i),B2(2,i)
 747         write (iout,*) 'CC'
 748         do j=1,2
 749           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 750         enddo
 751         write(iout,*) 'DD'
 752         do j=1,2
 753           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 754         enddo
 755         write(iout,*) 'EE'
 756         do j=1,2
 757           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 758         enddo
 759       enddo
 760       endif
 761 C
 762 C Read electrostatic-interaction parameters
 763 C
 764       if (lprint) then
 765         write (iout,*)
 766         write (iout,'(/a)') 'Electrostatic interaction constants:'
 767         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 768      &            'IT','JT','APP','BPP','AEL6','AEL3'
 769       endif
 770       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 771       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 772       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 773       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 774       close (ielep)
 775       do i=1,2
 776         do j=1,2
 777         rri=rpp(i,j)**6
 778         app (i,j)=epp(i,j)*rri*rri
 779         bpp (i,j)=-2.0D0*epp(i,j)*rri
 780         ael6(i,j)=elpp6(i,j)*4.2D0**6
 781         ael3(i,j)=elpp3(i,j)*4.2D0**3
 782         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 783      &                    ael6(i,j),ael3(i,j)
 784         enddo
 785       enddo
 786 C
 787 C Read side-chain interaction parameters.
 788 C
 789       read (isidep,*,end=117,err=117) ipot,expon
 790       if (ipot.lt.1 .or. ipot.gt.5) then
 791         write (iout,'(2a)') 'Error while reading SC interaction',
 792      &               'potential file - unknown potential type.'
 793 #ifdef MPI
 794         call MPI_Finalize(Ierror)
 795 #endif
 796         stop
 797       endif
 798       expon2=expon/2
 799       if(me.eq.king)
 800      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 801      & ', exponents are ',expon,2*expon
 802       goto (10,20,30,30,40) ipot
 803 C----------------------- LJ potential ---------------------------------
 804    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 805      &   (sigma0(i),i=1,ntyp)
 806       if (lprint) then
 807         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 808         write (iout,'(a/)') 'The epsilon array:'
 809         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 810         write (iout,'(/a)') 'One-body parameters:'
 811         write (iout,'(a,4x,a)') 'residue','sigma'
 812         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 813       endif
 814       goto 50
 815 C----------------------- LJK potential --------------------------------
 816    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 817      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 818       if (lprint) then
 819         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 820         write (iout,'(a/)') 'The epsilon array:'
 821         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 822         write (iout,'(/a)') 'One-body parameters:'
 823         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 824         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 825      &        i=1,ntyp)
 826       endif
 827       goto 50
 828 C---------------------- GB or BP potential -----------------------------
 829    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 830      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 831      &  (alp(i),i=1,ntyp)
 832 C For the GB potential convert sigma'**2 into chi'
 833       if (ipot.eq.4) then
 834         do i=1,ntyp
 835           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 836         enddo
 837       endif
 838       if (lprint) then
 839         write (iout,'(/a/)') 'Parameters of the BP potential:'
 840         write (iout,'(a/)') 'The epsilon array:'
 841         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 842         write (iout,'(/a)') 'One-body parameters:'
 843         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 844      &       '    chip  ','    alph  '
 845         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 846      &                     chip(i),alp(i),i=1,ntyp)
 847       endif
 848       goto 50
 849 C--------------------- GBV potential -----------------------------------
 850    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 851      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 852      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 853       if (lprint) then
 854         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 855         write (iout,'(a/)') 'The epsilon array:'
 856         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 857         write (iout,'(/a)') 'One-body parameters:'
 858         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 859      &      's||/s_|_^2','    chip  ','    alph  '
 860         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 861      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 862       endif
 863    50 continue
 864       close (isidep)
 865 C-----------------------------------------------------------------------
 866 C Calculate the "working" parameters of SC interactions.
 867       do i=2,ntyp
 868         do j=1,i-1
 869           eps(i,j)=eps(j,i)
 870         enddo
 871       enddo
 872       do i=1,ntyp
 873         do j=i,ntyp
 874           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 875           sigma(j,i)=sigma(i,j)
 876           rs0(i,j)=dwa16*sigma(i,j)
 877           rs0(j,i)=rs0(i,j)
 878         enddo
 879       enddo
 880       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 881      & 'Working parameters of the SC interactions:',
 882      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 883      & '  chi1   ','   chi2   '
 884       do i=1,ntyp
 885         do j=i,ntyp
 886           epsij=eps(i,j)
 887           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 888             rrij=sigma(i,j)
 889           else
 890             rrij=rr0(i)+rr0(j)
 891           endif
 892           r0(i,j)=rrij
 893           r0(j,i)=rrij
 894           rrij=rrij**expon
 895           epsij=eps(i,j)
 896           sigeps=dsign(1.0D0,epsij)
 897           epsij=dabs(epsij)
 898           aa(i,j)=epsij*rrij*rrij
 899           bb(i,j)=-sigeps*epsij*rrij
 900           aa(j,i)=aa(i,j)
 901           bb(j,i)=bb(i,j)
 902           if (ipot.gt.2) then
 903             sigt1sq=sigma0(i)**2
 904             sigt2sq=sigma0(j)**2
 905             sigii1=sigii(i)
 906             sigii2=sigii(j)
 907             ratsig1=sigt2sq/sigt1sq
 908             ratsig2=1.0D0/ratsig1
 909             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 910             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 911             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 912           else
 913             rsum_max=sigma(i,j)
 914           endif
 915 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 916             sigmaii(i,j)=rsum_max
 917             sigmaii(j,i)=rsum_max
 918 c         else
 919 c           sigmaii(i,j)=r0(i,j)
 920 c           sigmaii(j,i)=r0(i,j)
 921 c         endif
 922 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 923           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 924             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 925             augm(i,j)=epsij*r_augm**(2*expon)
 926 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 927             augm(j,i)=augm(i,j)
 928           else
 929             augm(i,j)=0.0D0
 930             augm(j,i)=0.0D0
 931           endif
 932           if (lprint) then
 933             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 934      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 935      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 936           endif
 937         enddo
 938       enddo
 939 #ifdef OLDSCP
 940 C
 941 C Define the SC-p interaction constants (hard-coded; old style)
 942 C
 943       do i=1,20
 944 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 945 C helix formation)
 946 c       aad(i,1)=0.3D0*4.0D0**12
 947 C Following line for constants currently implemented
 948 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 949         aad(i,1)=1.5D0*4.0D0**12
 950 c       aad(i,1)=0.17D0*5.6D0**12
 951         aad(i,2)=aad(i,1)
 952 C "Soft" SC-p repulsion
 953         bad(i,1)=0.0D0
 954 C Following line for constants currently implemented
 955 c       aad(i,1)=0.3D0*4.0D0**6
 956 C "Hard" SC-p repulsion
 957         bad(i,1)=3.0D0*4.0D0**6
 958 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 959         bad(i,2)=bad(i,1)
 960 c       aad(i,1)=0.0D0
 961 c       aad(i,2)=0.0D0
 962 c       bad(i,1)=1228.8D0
 963 c       bad(i,2)=1228.8D0
 964       enddo
 965 #else
 966 C
 967 C 8/9/01 Read the SC-p interaction constants from file
 968 C
 969       do i=1,ntyp
 970         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
 971       enddo
 972       do i=1,ntyp
 973         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 974         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 975         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 976         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 977       enddo
 978
 979       if (lprint) then
 980         write (iout,*) "Parameters of SC-p interactions:"
 981         do i=1,20
 982           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 983      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 984         enddo
 985       endif
 986 #endif
 987 C
 988 C Define the constants of the disulfide bridge
 989 C
 990       ebr=-5.50D0
 991 c
 992 c Old arbitrary potential - commented out.
 993 c
 994 c      dbr= 4.20D0
 995 c      fbr= 3.30D0
 996 c
 997 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 998 c energy surface of diethyl disulfide.
 999 c A. Liwo and U. Kozlowska, 11/24/03
1000 c
1001       D0CM = 3.78d0
1002       AKCM = 15.1d0
1003       AKTH = 11.0d0
1004       AKCT = 12.0d0
1005       V1SS =-1.08d0
1006       V2SS = 7.61d0
1007       V3SS = 13.7d0
1008 c      akcm=0.0d0
1009 c      akth=0.0d0
1010 c      akct=0.0d0
1011 c      v1ss=0.0d0
1012 c      v2ss=0.0d0
1013 c      v3ss=0.0d0
1014
1015       if(me.eq.king) then
1016       write (iout,'(/a)') "Disulfide bridge parameters:"
1017       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1018       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1019       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1020       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1021      &  ' v3ss:',v3ss
1022       endif
1023       return
1024   111 write (iout,*) "Error reading bending energy parameters."
1025       goto 999
1026   112 write (iout,*) "Error reading rotamer energy parameters."
1027       goto 999
1028   113 write (iout,*) "Error reading torsional energy parameters."
1029       goto 999
1030   114 write (iout,*) "Error reading double torsional energy parameters."
1031       goto 999
1032   115 write (iout,*)
1033      &  "Error reading cumulant (multibody energy) parameters."
1034       goto 999
1035   116 write (iout,*) "Error reading electrostatic energy parameters."
1036       goto 999
1037   117 write (iout,*) "Error reading side chain interaction parameters."
1038       goto 999
1039   118 write (iout,*) "Error reading SCp interaction parameters."
1040       goto 999
1041   119 write (iout,*) "Error reading SCCOR parameters"
1042   999 continue
1043 #ifdef MPI
1044       call MPI_Finalize(Ierror)
1045 #endif
1046       stop
1047       return
1048       end
1049
1050
1051       subroutine getenv_loc(var, val)
1052       character(*) var, val
1053
1054 #ifdef WINIFL
1055       character(2000) line
1056       external ilen
1057
1058       open (196,file='env',status='old',readonly,shared)
1059       iread=0
1060 c      write(*,*)'looking for ',var
1061 10    read(196,*,err=11,end=11)line
1062       iread=index(line,var)
1063 c      write(*,*)iread,' ',var,' ',line
1064       if (iread.eq.0) go to 10
1065 c      write(*,*)'---> ',line
1066 11    continue
1067       if(iread.eq.0) then
1068 c       write(*,*)'CHUJ'
1069        val=''
1070       else
1071        iread=iread+ilen(var)+1
1072        read (line(iread:),*,err=12,end=12) val
1073 c       write(*,*)'OK: ',var,' = ',val
1074       endif
1075       close(196)
1076       return
1077 12    val=''
1078       close(196)
1079 #elif (defined CRAY)
1080       integer lennam,lenval,ierror
1081 c
1082 c        getenv using a POSIX call, useful on the T3D
1083 c        Sept 1996, comment out error check on advice of H. Pritchard
1084 c
1085       lennam = len(var)
1086       if(lennam.le.0) stop '--error calling getenv--'
1087       call pxfgetenv(var,lennam,val,lenval,ierror)
1088 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1089 #else
1090       call getenv(var,val)
1091 #endif
1092
1093       return
1094       end