source/unres/src_MD/parmread.F

   1       subroutine parmread
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6 C Important! Energy-term weights ARE NOT read here; they are read from the
   7 C main input file instead, because NO defaults have yet been set for these
   8 C parameters.
   9 C
  10       implicit real*8 (a-h,o-z)
  11       include 'DIMENSIONS'
  12 #ifdef MPI
  13       include "mpif.h"
  14       integer IERROR
  15 #endif
  16       include 'COMMON.IOUNITS'
  17       include 'COMMON.CHAIN'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.GEO'
  20       include 'COMMON.LOCAL'
  21       include 'COMMON.TORSION'
  22       include 'COMMON.SCCOR'
  23       include 'COMMON.SCROT'
  24       include 'COMMON.FFIELD'
  25       include 'COMMON.NAMES'
  26       include 'COMMON.SBRIDGE'
  27       include 'COMMON.MD'
  28       include 'COMMON.SETUP'
  29       character*1 t1,t2,t3
  30       character*1 onelett(4) /"G","A","P","D"/
  31       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       enddo
 324       write (2,*) "End reading THETA_PDB"
 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+1,ntortyp-1
 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)=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 c              write (iout,*) "OK onelett",
 567 c     &         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 Martix 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 Martix 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      &           v1s(1,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 Read of Side-chain backbone correlation parameters
 656 C Modified 11 May 2012 by Adasko
 657 CCC
 658 C
 659       read (isccor,*,end=113,err=113) nsccortyp
 660 #ifdef SCCORPDB
 661       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
 662 #else
 663       read (isccor,*,end=113,err=113) (isccortyp(i),i=-ntyp,ntyp)
 664 #endif
 665 c      write (iout,*) 'ntortyp',ntortyp
 666       maxinter=3
 667 cc maxinter is maximum interaction sites
 668       do l=1,maxinter
 669       do i=1,nsccortyp
 670         do j=1,nsccortyp
 671           read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
 672           v0ijsccor=0.0d0
 673           si=-1.0d0
 674
 675           do k=1,nterm_sccor(i,j)
 676             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
 677      &    ,v2sccor(k,l,i,j)
 678             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 679             si=-si
 680           enddo
 681           do k=1,nlor_sccor(i,j)
 682             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 683      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 684             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 685      &(1+vlor3sccor(k,i,j)**2)
 686           enddo
 687           v0sccor(i,j)=v0ijsccor
 688         enddo
 689       enddo
 690       enddo
 691       close (isccor)
 692
 693       if (lprint) then
 694         write (iout,'(/a/)') 'Torsional constants:'
 695         do i=1,nsccortyp
 696           do j=1,nsccortyp
 697             write (iout,*) 'ityp',i,' jtyp',j
 698             write (iout,*) 'Fourier constants'
 699             do k=1,nterm_sccor(i,j)
 700       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 701             enddo
 702             write (iout,*) 'Lorenz constants'
 703             do k=1,nlor_sccor(i,j)
 704               write (iout,'(3(1pe15.5))')
 705      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 706             enddo
 707           enddo
 708         enddo
 709       endif
 710 C
 711 C
 712 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 713 C         interaction energy of the Gly, Ala, and Pro prototypes.
 714 C
 715       if (lprint) then
 716         write (iout,*)
 717         write (iout,*) "Coefficients of the cumulants"
 718       endif
 719       read (ifourier,*) nloctyp
 720       do i=0,nloctyp-1
 721         read (ifourier,*,end=115,err=115)
 722         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 723         if (lprint) then
 724         write (iout,*) 'Type',i
 725         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 726         endif
 727         B1(1,i)  = b(3)
 728         B1(2,i)  = b(5)
 729         B1(1,-i) = b(3)
 730         B1(2,-i) = -b(5)
 731 c        b1(1,i)=0.0d0
 732 c        b1(2,i)=0.0d0
 733         B1tilde(1,i) = b(3)
 734         B1tilde(2,i) =-b(5)
 735         B1tilde(1,-i) =-b(3)
 736         B1tilde(2,-i) =b(5)
 737 c        b1tilde(1,i)=0.0d0
 738 c        b1tilde(2,i)=0.0d0
 739         B2(1,i)  = b(2)
 740         B2(2,i)  = b(4)
 741         B2(1,-i)  =b(2)
 742         B2(2,-i)  =-b(4)
 743
 744 c        b2(1,i)=0.0d0
 745 c        b2(2,i)=0.0d0
 746         CC(1,1,i)= b(7)
 747         CC(2,2,i)=-b(7)
 748         CC(2,1,i)= b(9)
 749         CC(1,2,i)= b(9)
 750         CC(1,1,-i)= b(7)
 751         CC(2,2,-i)=-b(7)
 752         CC(2,1,-i)=-b(9)
 753         CC(1,2,-i)=-b(9)
 754 c        CC(1,1,i)=0.0d0
 755 c        CC(2,2,i)=0.0d0
 756 c        CC(2,1,i)=0.0d0
 757 c        CC(1,2,i)=0.0d0
 758         Ctilde(1,1,i)=b(7)
 759         Ctilde(1,2,i)=b(9)
 760         Ctilde(2,1,i)=-b(9)
 761         Ctilde(2,2,i)=b(7)
 762         Ctilde(1,1,-i)=b(7)
 763         Ctilde(1,2,-i)=-b(9)
 764         Ctilde(2,1,-i)=b(9)
 765         Ctilde(2,2,-i)=b(7)
 766
 767 c        Ctilde(1,1,i)=0.0d0
 768 c        Ctilde(1,2,i)=0.0d0
 769 c        Ctilde(2,1,i)=0.0d0
 770 c        Ctilde(2,2,i)=0.0d0
 771         DD(1,1,i)= b(6)
 772         DD(2,2,i)=-b(6)
 773         DD(2,1,i)= b(8)
 774         DD(1,2,i)= b(8)
 775         DD(1,1,-i)= b(6)
 776         DD(2,2,-i)=-b(6)
 777         DD(2,1,-i)=-b(8)
 778         DD(1,2,-i)=-b(8)
 779 c        DD(1,1,i)=0.0d0
 780 c        DD(2,2,i)=0.0d0
 781 c        DD(2,1,i)=0.0d0
 782 c        DD(1,2,i)=0.0d0
 783         Dtilde(1,1,i)=b(6)
 784         Dtilde(1,2,i)=b(8)
 785         Dtilde(2,1,i)=-b(8)
 786         Dtilde(2,2,i)=b(6)
 787         Dtilde(1,1,-i)=b(6)
 788         Dtilde(1,2,-i)=-b(8)
 789         Dtilde(2,1,-i)=b(8)
 790         Dtilde(2,2,-i)=b(6)
 791
 792 c        Dtilde(1,1,i)=0.0d0
 793 c        Dtilde(1,2,i)=0.0d0
 794 c        Dtilde(2,1,i)=0.0d0
 795 c        Dtilde(2,2,i)=0.0d0
 796         EE(1,1,i)= b(10)+b(11)
 797         EE(2,2,i)=-b(10)+b(11)
 798         EE(2,1,i)= b(12)-b(13)
 799         EE(1,2,i)= b(12)+b(13)
 800         EE(1,1,-i)= b(10)+b(11)
 801         EE(2,2,-i)=-b(10)+b(11)
 802         EE(2,1,-i)=-b(12)+b(13)
 803         EE(1,2,-i)=-b(12)-b(13)
 804
 805 c        ee(1,1,i)=1.0d0
 806 c        ee(2,2,i)=1.0d0
 807 c        ee(2,1,i)=0.0d0
 808 c        ee(1,2,i)=0.0d0
 809 c        ee(2,1,i)=ee(1,2,i)
 810       enddo
 811       if (lprint) then
 812       do i=1,nloctyp
 813         write (iout,*) 'Type',i
 814         write (iout,*) 'B1'
 815         write(iout,*) B1(1,i),B1(2,i)
 816         write (iout,*) 'B2'
 817         write(iout,*) B2(1,i),B2(2,i)
 818         write (iout,*) 'CC'
 819         do j=1,2
 820           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 821         enddo
 822         write(iout,*) 'DD'
 823         do j=1,2
 824           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 825         enddo
 826         write(iout,*) 'EE'
 827         do j=1,2
 828           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 829         enddo
 830       enddo
 831       endif
 832 C
 833 C Read electrostatic-interaction parameters
 834 C
 835       if (lprint) then
 836         write (iout,*)
 837         write (iout,'(/a)') 'Electrostatic interaction constants:'
 838         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 839      &            'IT','JT','APP','BPP','AEL6','AEL3'
 840       endif
 841       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 842       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 843       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 844       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 845       close (ielep)
 846       do i=1,2
 847         do j=1,2
 848         rri=rpp(i,j)**6
 849         app (i,j)=epp(i,j)*rri*rri
 850         bpp (i,j)=-2.0D0*epp(i,j)*rri
 851         ael6(i,j)=elpp6(i,j)*4.2D0**6
 852         ael3(i,j)=elpp3(i,j)*4.2D0**3
 853         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 854      &                    ael6(i,j),ael3(i,j)
 855         enddo
 856       enddo
 857 C
 858 C Read side-chain interaction parameters.
 859 C
 860       read (isidep,*,end=117,err=117) ipot,expon
 861       if (ipot.lt.1 .or. ipot.gt.5) then
 862         write (iout,'(2a)') 'Error while reading SC interaction',
 863      &               'potential file - unknown potential type.'
 864 #ifdef MPI
 865         call MPI_Finalize(Ierror)
 866 #endif
 867         stop
 868       endif
 869       expon2=expon/2
 870       if(me.eq.king)
 871      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 872      & ', exponents are ',expon,2*expon
 873       goto (10,20,30,30,40) ipot
 874 C----------------------- LJ potential ---------------------------------
 875    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 876      &   (sigma0(i),i=1,ntyp)
 877       if (lprint) then
 878         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 879         write (iout,'(a/)') 'The epsilon array:'
 880         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 881         write (iout,'(/a)') 'One-body parameters:'
 882         write (iout,'(a,4x,a)') 'residue','sigma'
 883         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 884       endif
 885       goto 50
 886 C----------------------- LJK potential --------------------------------
 887    20 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)
 889       if (lprint) then
 890         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 891         write (iout,'(a/)') 'The epsilon array:'
 892         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 893         write (iout,'(/a)') 'One-body parameters:'
 894         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 895         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 896      &        i=1,ntyp)
 897       endif
 898       goto 50
 899 C---------------------- GB or BP potential -----------------------------
 900    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 901      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 902      &  (alp(i),i=1,ntyp)
 903 C For the GB potential convert sigma'**2 into chi'
 904       if (ipot.eq.4) then
 905         do i=1,ntyp
 906           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 907         enddo
 908       endif
 909       if (lprint) then
 910         write (iout,'(/a/)') 'Parameters of the BP potential:'
 911         write (iout,'(a/)') 'The epsilon array:'
 912         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 913         write (iout,'(/a)') 'One-body parameters:'
 914         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 915      &       '    chip  ','    alph  '
 916         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 917      &                     chip(i),alp(i),i=1,ntyp)
 918       endif
 919       goto 50
 920 C--------------------- GBV potential -----------------------------------
 921    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 922      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 923      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 924       if (lprint) then
 925         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 926         write (iout,'(a/)') 'The epsilon array:'
 927         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 928         write (iout,'(/a)') 'One-body parameters:'
 929         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 930      &      's||/s_|_^2','    chip  ','    alph  '
 931         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 932      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 933       endif
 934    50 continue
 935       close (isidep)
 936 C-----------------------------------------------------------------------
 937 C Calculate the "working" parameters of SC interactions.
 938       do i=2,ntyp
 939         do j=1,i-1
 940           eps(i,j)=eps(j,i)
 941         enddo
 942       enddo
 943       do i=1,ntyp
 944         do j=i,ntyp
 945           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 946           sigma(j,i)=sigma(i,j)
 947           rs0(i,j)=dwa16*sigma(i,j)
 948           rs0(j,i)=rs0(i,j)
 949         enddo
 950       enddo
 951       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 952      & 'Working parameters of the SC interactions:',
 953      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 954      & '  chi1   ','   chi2   '
 955       do i=1,ntyp
 956         do j=i,ntyp
 957           epsij=eps(i,j)
 958           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 959             rrij=sigma(i,j)
 960           else
 961             rrij=rr0(i)+rr0(j)
 962           endif
 963           r0(i,j)=rrij
 964           r0(j,i)=rrij
 965           rrij=rrij**expon
 966           epsij=eps(i,j)
 967           sigeps=dsign(1.0D0,epsij)
 968           epsij=dabs(epsij)
 969           aa(i,j)=epsij*rrij*rrij
 970           bb(i,j)=-sigeps*epsij*rrij
 971           aa(j,i)=aa(i,j)
 972           bb(j,i)=bb(i,j)
 973           if (ipot.gt.2) then
 974             sigt1sq=sigma0(i)**2
 975             sigt2sq=sigma0(j)**2
 976             sigii1=sigii(i)
 977             sigii2=sigii(j)
 978             ratsig1=sigt2sq/sigt1sq
 979             ratsig2=1.0D0/ratsig1
 980             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 981             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 982             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 983           else
 984             rsum_max=sigma(i,j)
 985           endif
 986 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 987             sigmaii(i,j)=rsum_max
 988             sigmaii(j,i)=rsum_max
 989 c         else
 990 c           sigmaii(i,j)=r0(i,j)
 991 c           sigmaii(j,i)=r0(i,j)
 992 c         endif
 993 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 994           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 995             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 996             augm(i,j)=epsij*r_augm**(2*expon)
 997 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 998             augm(j,i)=augm(i,j)
 999           else
1000             augm(i,j)=0.0D0
1001             augm(j,i)=0.0D0
1002           endif
1003           if (lprint) then
1004             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1005      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1006      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1007           endif
1008         enddo
1009       enddo
1010 #ifdef OLDSCP
1011 C
1012 C Define the SC-p interaction constants (hard-coded; old style)
1013 C
1014       do i=1,20
1015 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1016 C helix formation)
1017 c       aad(i,1)=0.3D0*4.0D0**12
1018 C Following line for constants currently implemented
1019 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1020         aad(i,1)=1.5D0*4.0D0**12
1021 c       aad(i,1)=0.17D0*5.6D0**12
1022         aad(i,2)=aad(i,1)
1023 C "Soft" SC-p repulsion
1024         bad(i,1)=0.0D0
1025 C Following line for constants currently implemented
1026 c       aad(i,1)=0.3D0*4.0D0**6
1027 C "Hard" SC-p repulsion
1028         bad(i,1)=3.0D0*4.0D0**6
1029 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1030         bad(i,2)=bad(i,1)
1031 c       aad(i,1)=0.0D0
1032 c       aad(i,2)=0.0D0
1033 c       bad(i,1)=1228.8D0
1034 c       bad(i,2)=1228.8D0
1035       enddo
1036 #else
1037 C
1038 C 8/9/01 Read the SC-p interaction constants from file
1039 C
1040       do i=1,ntyp
1041         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1042       enddo
1043       do i=1,ntyp
1044         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1045         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1046         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1047         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1048       enddo
1049
1050       if (lprint) then
1051         write (iout,*) "Parameters of SC-p interactions:"
1052         do i=1,20
1053           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1054      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1055         enddo
1056       endif
1057 #endif
1058 C
1059 C Define the constants of the disulfide bridge
1060 C
1061       ebr=-5.50D0
1062 c
1063 c Old arbitrary potential - commented out.
1064 c
1065 c      dbr= 4.20D0
1066 c      fbr= 3.30D0
1067 c
1068 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1069 c energy surface of diethyl disulfide.
1070 c A. Liwo and U. Kozlowska, 11/24/03
1071 c
1072       D0CM = 3.78d0
1073       AKCM = 15.1d0
1074       AKTH = 11.0d0
1075       AKCT = 12.0d0
1076       V1SS =-1.08d0
1077       V2SS = 7.61d0
1078       V3SS = 13.7d0
1079 c      akcm=0.0d0
1080 c      akth=0.0d0
1081 c      akct=0.0d0
1082 c      v1ss=0.0d0
1083 c      v2ss=0.0d0
1084 c      v3ss=0.0d0
1085
1086       if(me.eq.king) then
1087       write (iout,'(/a)') "Disulfide bridge parameters:"
1088       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1089       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1090       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1091       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1092      &  ' v3ss:',v3ss
1093       endif
1094       return
1095   111 write (iout,*) "Error reading bending energy parameters."
1096       goto 999
1097   112 write (iout,*) "Error reading rotamer energy parameters."
1098       goto 999
1099   113 write (iout,*) "Error reading torsional energy parameters."
1100       goto 999
1101   114 write (iout,*) "Error reading double torsional energy parameters."
1102       goto 999
1103   115 write (iout,*)
1104      &  "Error reading cumulant (multibody energy) parameters."
1105       goto 999
1106   116 write (iout,*) "Error reading electrostatic energy parameters."
1107       goto 999
1108   117 write (iout,*) "Error reading side chain interaction parameters."
1109       goto 999
1110   118 write (iout,*) "Error reading SCp interaction parameters."
1111       goto 999
1112   119 write (iout,*) "Error reading SCCOR parameters"
1113   999 continue
1114 #ifdef MPI
1115       call MPI_Finalize(Ierror)
1116 #endif
1117       stop
1118       return
1119       end
1120
1121
1122       subroutine getenv_loc(var, val)
1123       character(*) var, val
1124
1125 #ifdef WINIFL
1126       character(2000) line
1127       external ilen
1128
1129       open (196,file='env',status='old',readonly,shared)
1130       iread=0
1131 c      write(*,*)'looking for ',var
1132 10    read(196,*,err=11,end=11)line
1133       iread=index(line,var)
1134 c      write(*,*)iread,' ',var,' ',line
1135       if (iread.eq.0) go to 10
1136 c      write(*,*)'---> ',line
1137 11    continue
1138       if(iread.eq.0) then
1139 c       write(*,*)'CHUJ'
1140        val=''
1141       else
1142        iread=iread+ilen(var)+1
1143        read (line(iread:),*,err=12,end=12) val
1144 c       write(*,*)'OK: ',var,' = ',val
1145       endif
1146       close(196)
1147       return
1148 12    val=''
1149       close(196)
1150 #elif (defined CRAY)
1151       integer lennam,lenval,ierror
1152 c
1153 c        getenv using a POSIX call, useful on the T3D
1154 c        Sept 1996, comment out error check on advice of H. Pritchard
1155 c
1156       lennam = len(var)
1157       if(lennam.le.0) stop '--error calling getenv--'
1158       call pxfgetenv(var,lennam,val,lenval,ierror)
1159 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1160 #else
1161       call getenv(var,val)
1162 #endif
1163
1164       return
1165       end