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