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