source/unres/src_MD/parmread.F

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