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