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,iblock)
 547               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 548      &         v2(k,i,j,iblock)
 549             enddo
 550             write (iout,*) 'Lorenz constants'
 551             do k=1,nlor(i,j,iblock)
 552               write (iout,'(3(1pe15.5))')
 553      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 554             enddo
 555           enddo
 556         enddo
 557       endif
 558 C
 559 C 6/23/01 Read parameters for double torsionals
 560 C
 561       do iblock=1,2
 562       do i=0,ntortyp-1
 563         do j=-ntortyp+1,ntortyp-1
 564           do k=-ntortyp+1,ntortyp-1
 565             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 566               write (iout,*) "OK onelett",
 567      &         i,j,k,t1,t2,t3
 568
 569             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 570      &        .or. t3.ne.toronelet(k)) then
 571               write (iout,*) "Error in double torsional parameter file",
 572      &         i,j,k,t1,t2,t3
 573 #ifdef MPI
 574               call MPI_Finalize(Ierror)
 575 #endif
 576                stop "Error in double torsional parameter file"
 577             endif
 578             read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
 579      &         ntermd_2(i,j,k,iblock)
 580             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 581             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 582             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
 583      &         ntermd_1(i,j,k,iblock))
 584             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
 585      &         ntermd_1(i,j,k,iblock))
 586             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
 587      &         ntermd_1(i,j,k,iblock))
 588             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
 589      &         ntermd_1(i,j,k,iblock))
 590 C Matrix of D parameters for one dimesional foureir series
 591             do l=1,ntermd_1(i,j,k,iblock)
 592              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 593              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 594              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 595              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 596 c            write(iout,*) "whcodze" ,
 597 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 598             enddo
 599             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
 600      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 601      &         v2s(m,l,i,j,k,iblock),
 602      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 603 C Matrix of D parameters for two dimesional fourier series
 604             do l=1,ntermd_2(i,j,k,iblock)
 605              do m=1,l-1
 606              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 607              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 608              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 609              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 610              enddo!m
 611             enddo!l
 612           enddo!k
 613         enddo!j
 614       enddo!i
 615       enddo!iblock
 616       if (lprint) then
 617       write (iout,*)
 618       write (iout,*) 'Constants for double torsionals'
 619       do iblock=1,2
 620       do i=1,ntortyp
 621         do j=-ntortyp,ntortyp
 622           do k=-ntortyp,ntortyp
 623             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 624      &        ' nsingle',ntermd_1(i,j,k,iblock),
 625      &        ' ndouble',ntermd_2(i,j,k,iblock)
 626             write (iout,*)
 627             write (iout,*) 'Single angles:'
 628             do l=1,ntermd_1(i,j,k,iblock)
 629               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 630      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 631      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
 632             enddo
 633             write (iout,*)
 634             write (iout,*) 'Pairs of angles:'
 635             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 636             do l=1,ntermd_2(i,j,k,iblock)
 637               write (iout,'(i5,20f10.5)')
 638      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 639             enddo
 640             write (iout,*)
 641             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 642             do l=1,ntermd_2(i,j,k,iblock)
 643               write (iout,'(i5,20f10.5)')
 644      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 645      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 646             enddo
 647             write (iout,*)
 648           enddo
 649         enddo
 650       enddo
 651       enddo
 652       endif
 653 #endif
 654 C
 655 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 656 C         correlation energies.
 657 C
 658       read (isccor,*,end=119,err=119) nterm_sccor
 659       do i=1,20
 660         do j=1,20
 661           read (isccor,'(a)')
 662           do k=1,nterm_sccor
 663             read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
 664      &        v2sccor(k,i,j)
 665           enddo
 666         enddo
 667       enddo
 668       close (isccor)
 669       if (lprint) then
 670         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 671         do i=1,20
 672           do j=1,20
 673             write (iout,*) 'ityp',i,' jtyp',j
 674             do k=1,nterm_sccor
 675               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
 676             enddo
 677           enddo
 678         enddo
 679       endif
 680 C
 681 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 682 C         interaction energy of the Gly, Ala, and Pro prototypes.
 683 C
 684       if (lprint) then
 685         write (iout,*)
 686         write (iout,*) "Coefficients of the cumulants"
 687       endif
 688       read (ifourier,*) nloctyp
 689       do i=1,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