source/maxlik/src-Fmatch/parmread.F.org

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