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