source/cluster/wham/src/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,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,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 Read of Side-chain backbone correlation parameters
 617 C Modified 11 May 2012 by Adasko
 618 CCC
 619 C
 620       read (isccor,*) nsccortyp
 621       read (isccor,*) (isccortyp(i),i=1,ntyp)
 622 c      write (iout,*) 'ntortyp',ntortyp
 623       do i=-ntyp,-1
 624         isccortyp(i)=-isccortyp(-i)
 625       enddo
 626       iscprol=isccortyp(20)
 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
 698       if (lprint) then
 699         write (iout,'(/a/)') 'Torsional constants:'
 700         do i=1,nsccortyp
 701           do j=1,nsccortyp
 702             write (iout,*) 'ityp',i,' jtyp',j
 703             write (iout,*) 'Fourier constants'
 704             do k=1,nterm_sccor(i,j)
 705               write (iout,'(2(1pe15.5))')
 706      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 707             enddo
 708             write (iout,*) 'Lorenz constants'
 709             do k=1,nlor_sccor(i,j)
 710               write (iout,'(3(1pe15.5))')
 711      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 712             enddo
 713           enddo
 714         enddo
 715       endif
 716
 717 C
 718 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 719 C         interaction energy of the Gly, Ala, and Pro prototypes.
 720 C
 721       read (ifourier,*) nloctyp
 722       do i=0,nloctyp-1
 723         read (ifourier,*)
 724         read (ifourier,*) (b(ii,i),ii=1,13)
 725         if (lprint) then
 726         write (iout,*) 'Type',i
 727         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 728         endif
 729         B1(1,i)  = b(3,i)
 730         B1(2,i)  = b(5,i)
 731         B1(1,-i) = b(3,i)
 732         B1(2,-i) = -b(5,i)
 733         B1tilde(1,i) = b(3,i)
 734         B1tilde(2,i) =-b(5,i)
 735         B1tilde(1,-i) =-b(3,i)
 736         B1tilde(2,-i) =b(5,i)
 737         B2(1,i)  = b(2,i)
 738         B2(2,i)  = b(4,i)
 739         B2(1,-i)  =b(2,i)
 740         B2(2,-i)  =-b(4,i)
 741         CC(1,1,i)= b(7,i)
 742         CC(2,2,i)=-b(7,i)
 743         CC(2,1,i)= b(9,i)
 744         CC(1,2,i)= b(9,i)
 745         CC(1,1,-i)= b(7,i)
 746         CC(2,2,-i)=-b(7,i)
 747         CC(2,1,-i)=-b(9,i)
 748         CC(1,2,-i)=-b(9,i)
 749         Ctilde(1,1,i)=b(7,i)
 750         Ctilde(1,2,i)=b(9,i)
 751         Ctilde(2,1,i)=-b(9,i)
 752         Ctilde(2,2,i)=b(7,i)
 753         Ctilde(1,1,-i)=b(7,i)
 754         Ctilde(1,2,-i)=-b(9,i)
 755         Ctilde(2,1,-i)=b(9,i)
 756         Ctilde(2,2,-i)=b(7,i)
 757         DD(1,1,i)= b(6,i)
 758         DD(2,2,i)=-b(6,i)
 759         DD(2,1,i)= b(8,i)
 760         DD(1,2,i)= b(8,i)
 761         DD(1,1,-i)= b(6,i)
 762         DD(2,2,-i)=-b(6,i)
 763         DD(2,1,-i)=-b(8,i)
 764         DD(1,2,-i)=-b(8,i)
 765         Dtilde(1,1,i)=b(6,i)
 766         Dtilde(1,2,i)=b(8,i)
 767         Dtilde(2,1,i)=-b(8,i)
 768         Dtilde(2,2,i)=b(6,i)
 769         Dtilde(1,1,-i)=b(6,i)
 770         Dtilde(1,2,-i)=-b(8,i)
 771         Dtilde(2,1,-i)=b(8,i)
 772         Dtilde(2,2,-i)=b(6,i)
 773         EE(1,1,i)= b(10,i)+b(11,i)
 774         EE(2,2,i)=-b(10,i)+b(11,i)
 775         EE(2,1,i)= b(12,i)-b(13,i)
 776         EE(1,2,i)= b(12,i)+b(13,i)
 777         EE(1,1,-i)= b(10,i)+b(11,i)
 778         EE(2,2,-i)=-b(10,i)+b(11,i)
 779         EE(2,1,-i)=-b(12,i)+b(13,i)
 780         EE(1,2,-i)=-b(12,i)-b(13,i)
 781       enddo
 782       if (lprint) then
 783       do i=1,nloctyp
 784         write (iout,*) 'Type',i
 785         write (iout,*) 'B1'
 786 c        write (iout,'(f10.5)') B1(:,i)
 787         write(iout,*) B1(1,i),B1(2,i)
 788         write (iout,*) 'B2'
 789 c        write (iout,'(f10.5)') B2(:,i)
 790         write(iout,*) B2(1,i),B2(2,i)
 791         write (iout,*) 'CC'
 792         do j=1,2
 793           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 794         enddo
 795         write(iout,*) 'DD'
 796         do j=1,2
 797           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 798         enddo
 799         write(iout,*) 'EE'
 800         do j=1,2
 801           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 802         enddo
 803       enddo
 804       endif
 805 C
 806 C Read electrostatic-interaction parameters
 807 C
 808       if (lprint) then
 809         write (iout,'(/a)') 'Electrostatic interaction constants:'
 810         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 811      &            'IT','JT','APP','BPP','AEL6','AEL3'
 812       endif
 813       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 814       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 815       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 816       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 817       close (ielep)
 818       do i=1,2
 819         do j=1,2
 820         rri=rpp(i,j)**6
 821         app (i,j)=epp(i,j)*rri*rri
 822         bpp (i,j)=-2.0D0*epp(i,j)*rri
 823         ael6(i,j)=elpp6(i,j)*4.2D0**6
 824         ael3(i,j)=elpp3(i,j)*4.2D0**3
 825         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 826      &                    ael6(i,j),ael3(i,j)
 827         enddo
 828       enddo
 829 C
 830 C Read side-chain interaction parameters.
 831 C
 832       read (isidep,*) ipot,expon
 833       if (ipot.lt.1 .or. ipot.gt.5) then
 834         write (iout,'(2a)') 'Error while reading SC interaction',
 835      &               'potential file - unknown potential type.'
 836         stop
 837       endif
 838       expon2=expon/2
 839       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 840      & ', exponents are ',expon,2*expon
 841       goto (10,20,30,30,40) ipot
 842 C----------------------- LJ potential ---------------------------------
 843    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 844       if (lprint) then
 845         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 846         write (iout,'(a/)') 'The epsilon array:'
 847         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 848         write (iout,'(/a)') 'One-body parameters:'
 849         write (iout,'(a,4x,a)') 'residue','sigma'
 850         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 851       endif
 852       goto 50
 853 C----------------------- LJK potential --------------------------------
 854    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 855      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 856       if (lprint) then
 857         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 858         write (iout,'(a/)') 'The epsilon array:'
 859         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 860         write (iout,'(/a)') 'One-body parameters:'
 861         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 862         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 863      &        i=1,ntyp)
 864       endif
 865       goto 50
 866 C---------------------- GB or BP potential -----------------------------
 867    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 868      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 869      &  (alp(i),i=1,ntyp)
 870 C For the GB potential convert sigma'**2 into chi'
 871       if (ipot.eq.4) then
 872         do i=1,ntyp
 873           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 874         enddo
 875       endif
 876       if (lprint) then
 877         write (iout,'(/a/)') 'Parameters of the BP potential:'
 878         write (iout,'(a/)') 'The epsilon array:'
 879         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 880         write (iout,'(/a)') 'One-body parameters:'
 881         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 882      &       '    chip  ','    alph  '
 883         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 884      &                     chip(i),alp(i),i=1,ntyp)
 885       endif
 886       goto 50
 887 C--------------------- GBV potential -----------------------------------
 888    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 889      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 890      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 891       if (lprint) then
 892         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 893         write (iout,'(a/)') 'The epsilon array:'
 894         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 895         write (iout,'(/a)') 'One-body parameters:'
 896         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 897      &      's||/s_|_^2','    chip  ','    alph  '
 898         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 899      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 900       endif
 901    50 continue
 902       close (isidep)
 903 C-----------------------------------------------------------------------
 904 C Calculate the "working" parameters of SC interactions.
 905       do i=2,ntyp
 906         do j=1,i-1
 907           eps(i,j)=eps(j,i)
 908         enddo
 909       enddo
 910       do i=1,ntyp
 911         do j=i,ntyp
 912           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 913           sigma(j,i)=sigma(i,j)
 914           rs0(i,j)=dwa16*sigma(i,j)
 915           rs0(j,i)=rs0(i,j)
 916         enddo
 917       enddo
 918       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 919      & 'Working parameters of the SC interactions:',
 920      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 921      & '  chi1   ','   chi2   '
 922       do i=1,ntyp
 923         do j=i,ntyp
 924           epsij=eps(i,j)
 925           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 926             rrij=sigma(i,j)
 927           else
 928             rrij=rr0(i)+rr0(j)
 929           endif
 930           r0(i,j)=rrij
 931           r0(j,i)=rrij
 932           rrij=rrij**expon
 933           epsij=eps(i,j)
 934           sigeps=dsign(1.0D0,epsij)
 935           epsij=dabs(epsij)
 936           aa(i,j)=epsij*rrij*rrij
 937           bb(i,j)=-sigeps*epsij*rrij
 938           aa(j,i)=aa(i,j)
 939           bb(j,i)=bb(i,j)
 940           if (ipot.gt.2) then
 941             sigt1sq=sigma0(i)**2
 942             sigt2sq=sigma0(j)**2
 943             sigii1=sigii(i)
 944             sigii2=sigii(j)
 945             ratsig1=sigt2sq/sigt1sq
 946             ratsig2=1.0D0/ratsig1
 947             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 948             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 949             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 950           else
 951             rsum_max=sigma(i,j)
 952           endif
 953 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 954             sigmaii(i,j)=rsum_max
 955             sigmaii(j,i)=rsum_max
 956 c         else
 957 c           sigmaii(i,j)=r0(i,j)
 958 c           sigmaii(j,i)=r0(i,j)
 959 c         endif
 960 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 961           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 962             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 963             augm(i,j)=epsij*r_augm**(2*expon)
 964 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 965             augm(j,i)=augm(i,j)
 966           else
 967             augm(i,j)=0.0D0
 968             augm(j,i)=0.0D0
 969           endif
 970           if (lprint) then
 971             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 972      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 973      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 974           endif
 975         enddo
 976       enddo
 977 C
 978 C Define the SC-p interaction constants
 979 C
 980       do i=1,20
 981         do j=1,2
 982           eps_scp(i,j)=-1.5d0
 983           rscp(i,j)=4.0d0
 984         enddo
 985       enddo
 986 #ifdef OLDSCP
 987       do i=1,20
 988 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 989 C helix formation)
 990 c       aad(i,1)=0.3D0*4.0D0**12
 991 C Following line for constants currently implemented
 992 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 993         aad(i,1)=1.5D0*4.0D0**12
 994 c       aad(i,1)=0.17D0*5.6D0**12
 995         aad(i,2)=aad(i,1)
 996 C "Soft" SC-p repulsion
 997         bad(i,1)=0.0D0
 998 C Following line for constants currently implemented
 999 c       aad(i,1)=0.3D0*4.0D0**6
1000 C "Hard" SC-p repulsion
1001         bad(i,1)=3.0D0*4.0D0**6
1002 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1003         bad(i,2)=bad(i,1)
1004 c       aad(i,1)=0.0D0
1005 c       aad(i,2)=0.0D0
1006 c       bad(i,1)=1228.8D0
1007 c       bad(i,2)=1228.8D0
1008       enddo
1009 #else
1010 C
1011 C 8/9/01 Read the SC-p interaction constants from file
1012 C
1013       do i=1,ntyp
1014         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1015       enddo
1016       do i=1,ntyp
1017         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1018         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1019         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1020         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1021       enddo
1022
1023       if (lprint) then
1024         write (iout,*) "Parameters of SC-p interactions:"
1025         do i=1,20
1026           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1027      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1028         enddo
1029       endif
1030 #endif
1031 C
1032 C Define the constants of the disulfide bridge
1033 C
1034 C      ebr=-5.50D0
1035 c
1036 c Old arbitrary potential - commented out.
1037 c
1038 c      dbr= 4.20D0
1039 c      fbr= 3.30D0
1040 c
1041 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1042 c energy surface of diethyl disulfide.
1043 c A. Liwo and U. Kozlowska, 11/24/03
1044 c
1045 C      D0CM = 3.78d0
1046 C      AKCM = 15.1d0
1047 C      AKTH = 11.0d0
1048 C      AKCT = 12.0d0
1049 C      V1SS =-1.08d0
1050 C      V2SS = 7.61d0
1051 C      V3SS = 13.7d0
1052       return
1053       end