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