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       if (lprint) then
 709       do i=1,nloctyp
 710         write (iout,*) 'Type',i
 711         write (iout,*) 'B1'
 712         write(iout,*) B1(1,i),B1(2,i)
 713         write (iout,*) 'B2'
 714         write(iout,*) B2(1,i),B2(2,i)
 715         write (iout,*) 'CC'
 716         do j=1,2
 717           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 718         enddo
 719         write(iout,*) 'DD'
 720         do j=1,2
 721           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 722         enddo
 723         write(iout,*) 'EE'
 724         do j=1,2
 725           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 726         enddo
 727       enddo
 728       endif
 729 C
 730 C Read electrostatic-interaction parameters
 731 C
 732       if (lprint) then
 733         write (iout,*)
 734         write (iout,'(/a)') 'Electrostatic interaction constants:'
 735         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 736      &            'IT','JT','APP','BPP','AEL6','AEL3'
 737       endif
 738       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 739       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 740       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 741       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 742       close (ielep)
 743       do i=1,2
 744         do j=1,2
 745         rri=rpp(i,j)**6
 746         app (i,j)=epp(i,j)*rri*rri
 747         bpp (i,j)=-2.0D0*epp(i,j)*rri
 748         ael6(i,j)=elpp6(i,j)*4.2D0**6
 749         ael3(i,j)=elpp3(i,j)*4.2D0**3
 750         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 751      &                    ael6(i,j),ael3(i,j)
 752         enddo
 753       enddo
 754 C
 755 C Read side-chain interaction parameters.
 756 C
 757       read (isidep,*,end=117,err=117) ipot,expon
 758       if (ipot.lt.1 .or. ipot.gt.5) then
 759         write (iout,'(2a)') 'Error while reading SC interaction',
 760      &               'potential file - unknown potential type.'
 761 #ifdef MPI
 762         call MPI_Finalize(Ierror)
 763 #endif
 764         stop
 765       endif
 766       expon2=expon/2
 767       if(me.eq.king)
 768      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 769      & ', exponents are ',expon,2*expon
 770       goto (10,20,30,30,40) ipot
 771 C----------------------- LJ potential ---------------------------------
 772    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 773      &   (sigma0(i),i=1,ntyp)
 774       if (lprint) then
 775         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 776         write (iout,'(a/)') 'The epsilon array:'
 777         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 778         write (iout,'(/a)') 'One-body parameters:'
 779         write (iout,'(a,4x,a)') 'residue','sigma'
 780         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 781       endif
 782       goto 50
 783 C----------------------- LJK potential --------------------------------
 784    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 785      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 786       if (lprint) then
 787         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 788         write (iout,'(a/)') 'The epsilon array:'
 789         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 790         write (iout,'(/a)') 'One-body parameters:'
 791         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 792         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 793      &        i=1,ntyp)
 794       endif
 795       goto 50
 796 C---------------------- GB or BP potential -----------------------------
 797    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 798      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 799      &  (alp(i),i=1,ntyp)
 800 C For the GB potential convert sigma'**2 into chi'
 801       if (ipot.eq.4) then
 802         do i=1,ntyp
 803           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 804         enddo
 805       endif
 806       if (lprint) then
 807         write (iout,'(/a/)') 'Parameters of the BP 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,4a)') 'residue','   sigma  ','s||/s_|_^2',
 812      &       '    chip  ','    alph  '
 813         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 814      &                     chip(i),alp(i),i=1,ntyp)
 815       endif
 816       goto 50
 817 C--------------------- GBV potential -----------------------------------
 818    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 819      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 820      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 821       if (lprint) then
 822         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 823         write (iout,'(a/)') 'The epsilon array:'
 824         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 825         write (iout,'(/a)') 'One-body parameters:'
 826         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 827      &      's||/s_|_^2','    chip  ','    alph  '
 828         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 829      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 830       endif
 831    50 continue
 832       close (isidep)
 833 C-----------------------------------------------------------------------
 834 C Calculate the "working" parameters of SC interactions.
 835       do i=2,ntyp
 836         do j=1,i-1
 837           eps(i,j)=eps(j,i)
 838         enddo
 839       enddo
 840       do i=1,ntyp
 841         do j=i,ntyp
 842           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 843           sigma(j,i)=sigma(i,j)
 844           rs0(i,j)=dwa16*sigma(i,j)
 845           rs0(j,i)=rs0(i,j)
 846         enddo
 847       enddo
 848       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 849      & 'Working parameters of the SC interactions:',
 850      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 851      & '  chi1   ','   chi2   '
 852       do i=1,ntyp
 853         do j=i,ntyp
 854           epsij=eps(i,j)
 855           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 856             rrij=sigma(i,j)
 857           else
 858             rrij=rr0(i)+rr0(j)
 859           endif
 860           r0(i,j)=rrij
 861           r0(j,i)=rrij
 862           rrij=rrij**expon
 863           epsij=eps(i,j)
 864           sigeps=dsign(1.0D0,epsij)
 865           epsij=dabs(epsij)
 866           aa(i,j)=epsij*rrij*rrij
 867           bb(i,j)=-sigeps*epsij*rrij
 868           aa(j,i)=aa(i,j)
 869           bb(j,i)=bb(i,j)
 870           if (ipot.gt.2) then
 871             sigt1sq=sigma0(i)**2
 872             sigt2sq=sigma0(j)**2
 873             sigii1=sigii(i)
 874             sigii2=sigii(j)
 875             ratsig1=sigt2sq/sigt1sq
 876             ratsig2=1.0D0/ratsig1
 877             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 878             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 879             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 880           else
 881             rsum_max=sigma(i,j)
 882           endif
 883 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 884             sigmaii(i,j)=rsum_max
 885             sigmaii(j,i)=rsum_max
 886 c         else
 887 c           sigmaii(i,j)=r0(i,j)
 888 c           sigmaii(j,i)=r0(i,j)
 889 c         endif
 890 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 891           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 892             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 893             augm(i,j)=epsij*r_augm**(2*expon)
 894 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 895             augm(j,i)=augm(i,j)
 896           else
 897             augm(i,j)=0.0D0
 898             augm(j,i)=0.0D0
 899           endif
 900           if (lprint) then
 901             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 902      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 903      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 904           endif
 905         enddo
 906       enddo
 907 #ifdef OLDSCP
 908 C
 909 C Define the SC-p interaction constants (hard-coded; old style)
 910 C
 911       do i=1,20
 912 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 913 C helix formation)
 914 c       aad(i,1)=0.3D0*4.0D0**12
 915 C Following line for constants currently implemented
 916 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 917         aad(i,1)=1.5D0*4.0D0**12
 918 c       aad(i,1)=0.17D0*5.6D0**12
 919         aad(i,2)=aad(i,1)
 920 C "Soft" SC-p repulsion
 921         bad(i,1)=0.0D0
 922 C Following line for constants currently implemented
 923 c       aad(i,1)=0.3D0*4.0D0**6
 924 C "Hard" SC-p repulsion
 925         bad(i,1)=3.0D0*4.0D0**6
 926 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 927         bad(i,2)=bad(i,1)
 928 c       aad(i,1)=0.0D0
 929 c       aad(i,2)=0.0D0
 930 c       bad(i,1)=1228.8D0
 931 c       bad(i,2)=1228.8D0
 932       enddo
 933 #else
 934 C
 935 C 8/9/01 Read the SC-p interaction constants from file
 936 C
 937       do i=1,ntyp
 938         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
 939       enddo
 940       do i=1,ntyp
 941         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 942         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 943         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 944         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 945       enddo
 946
 947       if (lprint) then
 948         write (iout,*) "Parameters of SC-p interactions:"
 949         do i=1,20
 950           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 951      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 952         enddo
 953       endif
 954 #endif
 955 C
 956 C Define the constants of the disulfide bridge
 957 C
 958       ebr=-5.50D0
 959 c
 960 c Old arbitrary potential - commented out.
 961 c
 962 c      dbr= 4.20D0
 963 c      fbr= 3.30D0
 964 c
 965 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 966 c energy surface of diethyl disulfide.
 967 c A. Liwo and U. Kozlowska, 11/24/03
 968 c
 969       D0CM = 3.78d0
 970       AKCM = 15.1d0
 971       AKTH = 11.0d0
 972       AKCT = 12.0d0
 973       V1SS =-1.08d0
 974       V2SS = 7.61d0
 975       V3SS = 13.7d0
 976 c      akcm=0.0d0
 977 c      akth=0.0d0
 978 c      akct=0.0d0
 979 c      v1ss=0.0d0
 980 c      v2ss=0.0d0
 981 c      v3ss=0.0d0
 982
 983       if(me.eq.king) then
 984       write (iout,'(/a)') "Disulfide bridge parameters:"
 985       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 986       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 987       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 988       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 989      &  ' v3ss:',v3ss
 990       endif
 991       return
 992   111 write (iout,*) "Error reading bending energy parameters."
 993       goto 999
 994   112 write (iout,*) "Error reading rotamer energy parameters."
 995       goto 999
 996   113 write (iout,*) "Error reading torsional energy parameters."
 997       goto 999
 998   114 write (iout,*) "Error reading double torsional energy parameters."
 999       goto 999
1000   115 write (iout,*)
1001      &  "Error reading cumulant (multibody energy) parameters."
1002       goto 999
1003   116 write (iout,*) "Error reading electrostatic energy parameters."
1004       goto 999
1005   117 write (iout,*) "Error reading side chain interaction parameters."
1006       goto 999
1007   118 write (iout,*) "Error reading SCp interaction parameters."
1008       goto 999
1009   119 write (iout,*) "Error reading SCCOR parameters"
1010   999 continue
1011 #ifdef MPI
1012       call MPI_Finalize(Ierror)
1013 #endif
1014       stop
1015       return
1016       end
1017
1018
1019       subroutine getenv_loc(var, val)
1020       character(*) var, val
1021
1022 #ifdef WINIFL
1023       character(2000) line
1024       external ilen
1025
1026       open (196,file='env',status='old',readonly,shared)
1027       iread=0
1028 c      write(*,*)'looking for ',var
1029 10    read(196,*,err=11,end=11)line
1030       iread=index(line,var)
1031 c      write(*,*)iread,' ',var,' ',line
1032       if (iread.eq.0) go to 10
1033 c      write(*,*)'---> ',line
1034 11    continue
1035       if(iread.eq.0) then
1036 c       write(*,*)'CHUJ'
1037        val=''
1038       else
1039        iread=iread+ilen(var)+1
1040        read (line(iread:),*,err=12,end=12) val
1041 c       write(*,*)'OK: ',var,' = ',val
1042       endif
1043       close(196)
1044       return
1045 12    val=''
1046       close(196)
1047 #elif (defined CRAY)
1048       integer lennam,lenval,ierror
1049 c
1050 c        getenv using a POSIX call, useful on the T3D
1051 c        Sept 1996, comment out error check on advice of H. Pritchard
1052 c
1053       lennam = len(var)
1054       if(lennam.le.0) stop '--error calling getenv--'
1055       call pxfgetenv(var,lennam,val,lenval,ierror)
1056 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1057 #else
1058       call getenv(var,val)
1059 #endif
1060
1061       return
1062       end