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