source/unres/src_CSA_DiL/parmread.F

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