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