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