source/unres/src_MD-M/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 #endif
 280       close(ithep)
 281 #ifdef CRYST_SC
 282 C
 283 C Read the parameters of the probability distribution/energy expression
 284 C of the side chains.
 285 C
 286       do i=1,ntyp
 287         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 288         if (i.eq.10) then
 289           dsc_inv(i)=0.0D0
 290         else
 291           dsc_inv(i)=1.0D0/dsc(i)
 292         endif
 293         if (i.ne.10) then
 294         do j=1,nlob(i)
 295           do k=1,3
 296             do l=1,3
 297               blower(l,k,j)=0.0D0
 298             enddo
 299           enddo
 300         enddo
 301         bsc(1,i)=0.0D0
 302         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 303      &    ((blower(k,l,1),l=1,k),k=1,3)
 304         do j=2,nlob(i)
 305           read (irotam,*,end=112,err=112) bsc(j,i)
 306           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 307      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 308         enddo
 309         do j=1,nlob(i)
 310           do k=1,3
 311             do l=1,k
 312               akl=0.0D0
 313               do m=1,3
 314                 akl=akl+blower(k,m,j)*blower(l,m,j)
 315               enddo
 316               gaussc(k,l,j,i)=akl
 317               gaussc(l,k,j,i)=akl
 318             enddo
 319           enddo
 320         enddo
 321         endif
 322       enddo
 323       close (irotam)
 324       if (lprint) then
 325         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 326         do i=1,ntyp
 327           nlobi=nlob(i)
 328           if (nlobi.gt.0) then
 329             if (LaTeX) then
 330               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 331      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 332                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 333      &                             'log h',(bsc(j,i),j=1,nlobi)
 334                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 335      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 336               do k=1,3
 337                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 338      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 339               enddo
 340             else
 341               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 342               write (iout,'(a,f10.4,4(16x,f10.4))')
 343      &                             'Center  ',(bsc(j,i),j=1,nlobi)
 344               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
 345      &           j=1,nlobi)
 346               write (iout,'(a)')
 347             endif
 348           endif
 349         enddo
 350       endif
 351 #else
 352 C
 353 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 354 C added by Urszula Kozlowska 07/11/2007
 355 C
 356       do i=1,ntyp
 357         read (irotam,*,end=112,err=112)
 358        if (i.eq.10) then
 359          read (irotam,*,end=112,err=112)
 360        else
 361          do j=1,65
 362            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 363          enddo
 364        endif
 365       enddo
 366 #endif
 367       close(irotam)
 368
 369 #ifdef CRYST_TOR
 370 C
 371 C Read torsional parameters in old format
 372 C
 373       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 374       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 375       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 376       do i=1,ntortyp
 377         do j=1,ntortyp
 378           read (itorp,'(a)')
 379           do k=1,nterm_old
 380             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 381           enddo
 382         enddo
 383       enddo
 384       close (itorp)
 385       if (lprint) then
 386         write (iout,'(/a/)') 'Torsional constants:'
 387         do i=1,ntortyp
 388           do j=1,ntortyp
 389             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 390             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 391           enddo
 392         enddo
 393       endif
 394 #else
 395 C
 396 C Read torsional parameters
 397 C
 398       read (itorp,*,end=113,err=113) ntortyp
 399       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 400       itortyp(ntyp1)=ntortyp+1
 401 c      write (iout,*) 'ntortyp',ntortyp
 402       do i=1,ntortyp
 403         do j=1,ntortyp
 404           read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
 405           v0ij=0.0d0
 406           si=-1.0d0
 407           do k=1,nterm(i,j)
 408             read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
 409             v0ij=v0ij+si*v1(k,i,j)
 410             si=-si
 411           enddo
 412           do k=1,nlor(i,j)
 413             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 414      &        vlor2(k,i,j),vlor3(k,i,j)
 415             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 416           enddo
 417           v0(i,j)=v0ij
 418         enddo
 419       enddo
 420       close (itorp)
 421       if (lprint) then
 422         write (iout,'(/a/)') 'Torsional constants:'
 423         do i=1,ntortyp
 424           do j=1,ntortyp
 425             write (iout,*) 'ityp',i,' jtyp',j
 426             write (iout,*) 'Fourier constants'
 427             do k=1,nterm(i,j)
 428               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 429             enddo
 430             write (iout,*) 'Lorenz constants'
 431             do k=1,nlor(i,j)
 432               write (iout,'(3(1pe15.5))')
 433      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 434             enddo
 435           enddo
 436         enddo
 437       endif
 438 C
 439 C 6/23/01 Read parameters for double torsionals
 440 C
 441       do i=1,ntortyp
 442         do j=1,ntortyp
 443           do k=1,ntortyp
 444             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 445             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 446      &        .or. t3.ne.onelett(k)) then
 447               write (iout,*) "Error in double torsional parameter file",
 448      &         i,j,k,t1,t2,t3
 449 #ifdef MPI
 450               call MPI_Finalize(Ierror)
 451 #endif
 452                stop "Error in double torsional parameter file"
 453             endif
 454             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
 455      &         ntermd_2(i,j,k)
 456             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
 457      &         ntermd_1(i,j,k))
 458             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
 459      &         ntermd_1(i,j,k))
 460             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
 461      &         ntermd_1(i,j,k))
 462             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
 463      &         ntermd_1(i,j,k))
 464             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
 465      &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
 466      &         m=1,l-1),l=1,ntermd_2(i,j,k))
 467           enddo
 468         enddo
 469       enddo
 470       if (lprint) then
 471       write (iout,*)
 472       write (iout,*) 'Constants for double torsionals'
 473       do i=1,ntortyp
 474         do j=1,ntortyp
 475           do k=1,ntortyp
 476             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 477      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 478             write (iout,*)
 479             write (iout,*) 'Single angles:'
 480             do l=1,ntermd_1(i,j,k)
 481               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 482      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 483      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 484             enddo
 485             write (iout,*)
 486             write (iout,*) 'Pairs of angles:'
 487             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 488             do l=1,ntermd_2(i,j,k)
 489               write (iout,'(i5,20f10.5)')
 490      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 491             enddo
 492             write (iout,*)
 493             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 494             do l=1,ntermd_2(i,j,k)
 495               write (iout,'(i5,20f10.5)')
 496      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 497             enddo
 498             write (iout,*)
 499           enddo
 500         enddo
 501       enddo
 502       endif
 503 #endif
 504 C Read of Side-chain backbone correlation parameters
 505 C Modified 11 May 2012 by Adasko
 506 CCC
 507 C
 508 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 509 C
 510       read (isccor,*,end=119,err=119) nsccortyp
 511 #ifdef SCCORPDB
 512       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 513       do i=-ntyp,-1
 514         isccortyp(i)=-isccortyp(-i)
 515       enddo
 516       iscprol=isccortyp(20)
 517 c      write (iout,*) 'ntortyp',ntortyp
 518       maxinter=3
 519 cc maxinter is maximum interaction sites
 520       do l=1,maxinter
 521       do i=1,nsccortyp
 522         do j=1,nsccortyp
 523           read (isccor,*,end=119,err=119)
 524      &nterm_sccor(i,j),nlor_sccor(i,j)
 525           v0ijsccor=0.0d0
 526           v0ijsccor1=0.0d0
 527           v0ijsccor2=0.0d0
 528           v0ijsccor3=0.0d0
 529           si=-1.0d0
 530           nterm_sccor(-i,j)=nterm_sccor(i,j)
 531           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 532           nterm_sccor(i,-j)=nterm_sccor(i,j)
 533           do k=1,nterm_sccor(i,j)
 534             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 535      &    ,v2sccor(k,l,i,j)
 536             if (j.eq.iscprol) then
 537              if (i.eq.isccortyp(10)) then
 538              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 539              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 540              else
 541              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 542      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 543              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 544      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 545              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 546              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 547              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 548              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 549              endif
 550             else
 551              if (i.eq.isccortyp(10)) then
 552              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 553              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 554              else
 555                if (j.eq.isccortyp(10)) then
 556              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 557              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 558                else
 559              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 560              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 561              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 562              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 563              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 564              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 565                 endif
 566                endif
 567             endif
 568             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 569             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 570             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 571             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 572             si=-si
 573           enddo
 574           do k=1,nlor_sccor(i,j)
 575             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 576      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 577             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 578      &(1+vlor3sccor(k,i,j)**2)
 579           enddo
 580           v0sccor(l,i,j)=v0ijsccor
 581           v0sccor(l,-i,j)=v0ijsccor1
 582           v0sccor(l,i,-j)=v0ijsccor2
 583           v0sccor(l,-i,-j)=v0ijsccor3
 584         enddo
 585       enddo
 586       enddo
 587       close (isccor)
 588 #else
 589       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
 590 c      write (iout,*) 'ntortyp',ntortyp
 591       maxinter=3
 592 cc maxinter is maximum interaction sites
 593       do l=1,maxinter
 594       do i=1,nsccortyp
 595         do j=1,nsccortyp
 596           read (isccor,*,end=113,err=113)
 597      & nterm_sccor(i,j),nlor_sccor(i,j)
 598           v0ijsccor=0.0d0
 599           si=-1.0d0
 600
 601           do k=1,nterm_sccor(i,j)
 602             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
 603      &    ,v2sccor(k,l,i,j)
 604             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 605             si=-si
 606           enddo
 607           do k=1,nlor_sccor(i,j)
 608             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 609      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 610             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 611      &(1+vlor3sccor(k,i,j)**2)
 612           enddo
 613           v0sccor(l,i,j)=v0ijsccor
 614         enddo
 615       enddo
 616       enddo
 617       close (isccor)
 618
 619 #endif
 620       if (lprint) then
 621         write (iout,'(/a/)') 'Torsional constants:'
 622         do i=1,nsccortyp
 623           do j=1,nsccortyp
 624             write (iout,*) 'ityp',i,' jtyp',j
 625             write (iout,*) 'Fourier constants'
 626             do k=1,nterm_sccor(i,j)
 627       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 628             enddo
 629             write (iout,*) 'Lorenz constants'
 630             do k=1,nlor_sccor(i,j)
 631               write (iout,'(3(1pe15.5))')
 632      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 633             enddo
 634           enddo
 635         enddo
 636       endif
 637
 638 C
 639 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 640 C         interaction energy of the Gly, Ala, and Pro prototypes.
 641 C
 642       if (lprint) then
 643         write (iout,*)
 644         write (iout,*) "Coefficients of the cumulants"
 645       endif
 646       read (ifourier,*) nloctyp
 647       do i=1,nloctyp
 648         read (ifourier,*,end=115,err=115)
 649         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 650         if (lprint) then
 651         write (iout,*) 'Type',i
 652         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 653         endif
 654         B1(1,i)  = b(3)
 655         B1(2,i)  = b(5)
 656 c        b1(1,i)=0.0d0
 657 c        b1(2,i)=0.0d0
 658         B1tilde(1,i) = b(3)
 659         B1tilde(2,i) =-b(5)
 660 c        b1tilde(1,i)=0.0d0
 661 c        b1tilde(2,i)=0.0d0
 662         B2(1,i)  = b(2)
 663         B2(2,i)  = b(4)
 664 c        b2(1,i)=0.0d0
 665 c        b2(2,i)=0.0d0
 666         CC(1,1,i)= b(7)
 667         CC(2,2,i)=-b(7)
 668         CC(2,1,i)= b(9)
 669         CC(1,2,i)= b(9)
 670 c        CC(1,1,i)=0.0d0
 671 c        CC(2,2,i)=0.0d0
 672 c        CC(2,1,i)=0.0d0
 673 c        CC(1,2,i)=0.0d0
 674         Ctilde(1,1,i)=b(7)
 675         Ctilde(1,2,i)=b(9)
 676         Ctilde(2,1,i)=-b(9)
 677         Ctilde(2,2,i)=b(7)
 678 c        Ctilde(1,1,i)=0.0d0
 679 c        Ctilde(1,2,i)=0.0d0
 680 c        Ctilde(2,1,i)=0.0d0
 681 c        Ctilde(2,2,i)=0.0d0
 682         DD(1,1,i)= b(6)
 683         DD(2,2,i)=-b(6)
 684         DD(2,1,i)= b(8)
 685         DD(1,2,i)= b(8)
 686 c        DD(1,1,i)=0.0d0
 687 c        DD(2,2,i)=0.0d0
 688 c        DD(2,1,i)=0.0d0
 689 c        DD(1,2,i)=0.0d0
 690         Dtilde(1,1,i)=b(6)
 691         Dtilde(1,2,i)=b(8)
 692         Dtilde(2,1,i)=-b(8)
 693         Dtilde(2,2,i)=b(6)
 694 c        Dtilde(1,1,i)=0.0d0
 695 c        Dtilde(1,2,i)=0.0d0
 696 c        Dtilde(2,1,i)=0.0d0
 697 c        Dtilde(2,2,i)=0.0d0
 698         EE(1,1,i)= b(10)+b(11)
 699         EE(2,2,i)=-b(10)+b(11)
 700         EE(2,1,i)= b(12)-b(13)
 701         EE(1,2,i)= b(12)+b(13)
 702 c        ee(1,1,i)=1.0d0
 703 c        ee(2,2,i)=1.0d0
 704 c        ee(2,1,i)=0.0d0
 705 c        ee(1,2,i)=0.0d0
 706 c        ee(2,1,i)=ee(1,2,i)
 707       enddo
 708 C      write(iout,*) "parm", B1(1,nloctyp+1),B1(2,nloctyp+1)
 709       if (lprint) then
 710       do i=1,nloctyp+1
 711         write (iout,*) 'Type',i
 712         write (iout,*) 'B1'
 713         write(iout,*) B1(1,i),B1(2,i)
 714         write (iout,*) 'B2'
 715         write(iout,*) B2(1,i),B2(2,i)
 716         write (iout,*) 'CC'
 717         do j=1,2
 718           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 719         enddo
 720         write(iout,*) 'DD'
 721         do j=1,2
 722           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 723         enddo
 724         write(iout,*) 'EE'
 725         do j=1,2
 726           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 727         enddo
 728       enddo
 729       endif
 730 C
 731 C Read electrostatic-interaction parameters
 732 C
 733       if (lprint) then
 734         write (iout,*)
 735         write (iout,'(/a)') 'Electrostatic interaction constants:'
 736         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 737      &            'IT','JT','APP','BPP','AEL6','AEL3'
 738       endif
 739       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 740       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 741       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 742       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 743       close (ielep)
 744       do i=1,2
 745         do j=1,2
 746         rri=rpp(i,j)**6
 747         app (i,j)=epp(i,j)*rri*rri
 748         bpp (i,j)=-2.0D0*epp(i,j)*rri
 749         ael6(i,j)=elpp6(i,j)*4.2D0**6
 750         ael3(i,j)=elpp3(i,j)*4.2D0**3
 751         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 752      &                    ael6(i,j),ael3(i,j)
 753         enddo
 754       enddo
 755 C
 756 C Read side-chain interaction parameters.
 757 C
 758       read (isidep,*,end=117,err=117) ipot,expon
 759       if (ipot.lt.1 .or. ipot.gt.5) then
 760         write (iout,'(2a)') 'Error while reading SC interaction',
 761      &               'potential file - unknown potential type.'
 762 #ifdef MPI
 763         call MPI_Finalize(Ierror)
 764 #endif
 765         stop
 766       endif
 767       expon2=expon/2
 768       if(me.eq.king)
 769      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 770      & ', exponents are ',expon,2*expon
 771       goto (10,20,30,30,40) ipot
 772 C----------------------- LJ potential ---------------------------------
 773    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 774      &   (sigma0(i),i=1,ntyp)
 775       if (lprint) then
 776         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 777         write (iout,'(a/)') 'The epsilon array:'
 778         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 779         write (iout,'(/a)') 'One-body parameters:'
 780         write (iout,'(a,4x,a)') 'residue','sigma'
 781         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 782       endif
 783       goto 50
 784 C----------------------- LJK potential --------------------------------
 785    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 786      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 787       if (lprint) then
 788         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 789         write (iout,'(a/)') 'The epsilon array:'
 790         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 791         write (iout,'(/a)') 'One-body parameters:'
 792         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 793         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 794      &        i=1,ntyp)
 795       endif
 796       goto 50
 797 C---------------------- GB or BP potential -----------------------------
 798    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 799      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 800      &  (alp(i),i=1,ntyp)
 801 C For the GB potential convert sigma'**2 into chi'
 802       if (ipot.eq.4) then
 803         do i=1,ntyp
 804           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 805         enddo
 806       endif
 807       if (lprint) then
 808         write (iout,'(/a/)') 'Parameters of the BP potential:'
 809         write (iout,'(a/)') 'The epsilon array:'
 810         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 811         write (iout,'(/a)') 'One-body parameters:'
 812         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 813      &       '    chip  ','    alph  '
 814         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 815      &                     chip(i),alp(i),i=1,ntyp)
 816       endif
 817       goto 50
 818 C--------------------- GBV potential -----------------------------------
 819    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 820      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 821      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 822       if (lprint) then
 823         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 824         write (iout,'(a/)') 'The epsilon array:'
 825         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 826         write (iout,'(/a)') 'One-body parameters:'
 827         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 828      &      's||/s_|_^2','    chip  ','    alph  '
 829         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 830      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 831       endif
 832    50 continue
 833       close (isidep)
 834 C-----------------------------------------------------------------------
 835 C Calculate the "working" parameters of SC interactions.
 836       do i=2,ntyp
 837         do j=1,i-1
 838           eps(i,j)=eps(j,i)
 839         enddo
 840       enddo
 841       do i=1,ntyp
 842         do j=i,ntyp
 843           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 844           sigma(j,i)=sigma(i,j)
 845           rs0(i,j)=dwa16*sigma(i,j)
 846           rs0(j,i)=rs0(i,j)
 847         enddo
 848       enddo
 849       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 850      & 'Working parameters of the SC interactions:',
 851      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 852      & '  chi1   ','   chi2   '
 853       do i=1,ntyp
 854         do j=i,ntyp
 855           epsij=eps(i,j)
 856           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 857             rrij=sigma(i,j)
 858           else
 859             rrij=rr0(i)+rr0(j)
 860           endif
 861           r0(i,j)=rrij
 862           r0(j,i)=rrij
 863           rrij=rrij**expon
 864           epsij=eps(i,j)
 865           sigeps=dsign(1.0D0,epsij)
 866           epsij=dabs(epsij)
 867           aa(i,j)=epsij*rrij*rrij
 868           bb(i,j)=-sigeps*epsij*rrij
 869           aa(j,i)=aa(i,j)
 870           bb(j,i)=bb(i,j)
 871           if (ipot.gt.2) then
 872             sigt1sq=sigma0(i)**2
 873             sigt2sq=sigma0(j)**2
 874             sigii1=sigii(i)
 875             sigii2=sigii(j)
 876             ratsig1=sigt2sq/sigt1sq
 877             ratsig2=1.0D0/ratsig1
 878             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 879             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 880             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 881           else
 882             rsum_max=sigma(i,j)
 883           endif
 884 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 885             sigmaii(i,j)=rsum_max
 886             sigmaii(j,i)=rsum_max
 887 c         else
 888 c           sigmaii(i,j)=r0(i,j)
 889 c           sigmaii(j,i)=r0(i,j)
 890 c         endif
 891 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 892           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 893             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 894             augm(i,j)=epsij*r_augm**(2*expon)
 895 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 896             augm(j,i)=augm(i,j)
 897           else
 898             augm(i,j)=0.0D0
 899             augm(j,i)=0.0D0
 900           endif
 901           if (lprint) then
 902             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 903      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 904      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 905           endif
 906         enddo
 907       enddo
 908 #ifdef OLDSCP
 909 C
 910 C Define the SC-p interaction constants (hard-coded; old style)
 911 C
 912       do i=1,20
 913 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 914 C helix formation)
 915 c       aad(i,1)=0.3D0*4.0D0**12
 916 C Following line for constants currently implemented
 917 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 918         aad(i,1)=1.5D0*4.0D0**12
 919 c       aad(i,1)=0.17D0*5.6D0**12
 920         aad(i,2)=aad(i,1)
 921 C "Soft" SC-p repulsion
 922         bad(i,1)=0.0D0
 923 C Following line for constants currently implemented
 924 c       aad(i,1)=0.3D0*4.0D0**6
 925 C "Hard" SC-p repulsion
 926         bad(i,1)=3.0D0*4.0D0**6
 927 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 928         bad(i,2)=bad(i,1)
 929 c       aad(i,1)=0.0D0
 930 c       aad(i,2)=0.0D0
 931 c       bad(i,1)=1228.8D0
 932 c       bad(i,2)=1228.8D0
 933       enddo
 934 #else
 935 C
 936 C 8/9/01 Read the SC-p interaction constants from file
 937 C
 938       do i=1,ntyp
 939         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
 940       enddo
 941       do i=1,ntyp
 942         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 943         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 944         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 945         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 946       enddo
 947
 948       if (lprint) then
 949         write (iout,*) "Parameters of SC-p interactions:"
 950         do i=1,20
 951           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 952      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 953         enddo
 954       endif
 955 #endif
 956 C
 957 C Define the constants of the disulfide bridge
 958 C
 959       ebr=-5.50D0
 960 c
 961 c Old arbitrary potential - commented out.
 962 c
 963 c      dbr= 4.20D0
 964 c      fbr= 3.30D0
 965 c
 966 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 967 c energy surface of diethyl disulfide.
 968 c A. Liwo and U. Kozlowska, 11/24/03
 969 c
 970       D0CM = 3.78d0
 971       AKCM = 15.1d0
 972       AKTH = 11.0d0
 973       AKCT = 12.0d0
 974       V1SS =-1.08d0
 975       V2SS = 7.61d0
 976       V3SS = 13.7d0
 977 c      akcm=0.0d0
 978 c      akth=0.0d0
 979 c      akct=0.0d0
 980 c      v1ss=0.0d0
 981 c      v2ss=0.0d0
 982 c      v3ss=0.0d0
 983
 984       if(me.eq.king) then
 985       write (iout,'(/a)') "Disulfide bridge parameters:"
 986       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 987       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 988       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 989       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 990      &  ' v3ss:',v3ss
 991       endif
 992       return
 993   111 write (iout,*) "Error reading bending energy parameters."
 994       goto 999
 995   112 write (iout,*) "Error reading rotamer energy parameters."
 996       goto 999
 997   113 write (iout,*) "Error reading torsional energy parameters."
 998       goto 999
 999   114 write (iout,*) "Error reading double torsional energy parameters."
1000       goto 999
1001   115 write (iout,*)
1002      &  "Error reading cumulant (multibody energy) parameters."
1003       goto 999
1004   116 write (iout,*) "Error reading electrostatic energy parameters."
1005       goto 999
1006   117 write (iout,*) "Error reading side chain interaction parameters."
1007       goto 999
1008   118 write (iout,*) "Error reading SCp interaction parameters."
1009       goto 999
1010   119 write (iout,*) "Error reading SCCOR parameters"
1011   999 continue
1012 #ifdef MPI
1013       call MPI_Finalize(Ierror)
1014 #endif
1015       stop
1016       return
1017       end
1018
1019
1020       subroutine getenv_loc(var, val)
1021       character(*) var, val
1022
1023 #ifdef WINIFL
1024       character(2000) line
1025       external ilen
1026
1027       open (196,file='env',status='old',readonly,shared)
1028       iread=0
1029 c      write(*,*)'looking for ',var
1030 10    read(196,*,err=11,end=11)line
1031       iread=index(line,var)
1032 c      write(*,*)iread,' ',var,' ',line
1033       if (iread.eq.0) go to 10
1034 c      write(*,*)'---> ',line
1035 11    continue
1036       if(iread.eq.0) then
1037 c       write(*,*)'CHUJ'
1038        val=''
1039       else
1040        iread=iread+ilen(var)+1
1041        read (line(iread:),*,err=12,end=12) val
1042 c       write(*,*)'OK: ',var,' = ',val
1043       endif
1044       close(196)
1045       return
1046 12    val=''
1047       close(196)
1048 #elif (defined CRAY)
1049       integer lennam,lenval,ierror
1050 c
1051 c        getenv using a POSIX call, useful on the T3D
1052 c        Sept 1996, comment out error check on advice of H. Pritchard
1053 c
1054       lennam = len(var)
1055       if(lennam.le.0) stop '--error calling getenv--'
1056       call pxfgetenv(var,lennam,val,lenval,ierror)
1057 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1058 #else
1059       call getenv(var,val)
1060 #endif
1061
1062       return
1063       end