source/wham/src/parmread.F

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