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       logical lprint
  27       dimension blower(3,3,maxlob)
  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       integer ilen
  33       external ilen
  34       character*16 key
  35       integer iparm
  36       double precision ip,mp
  37 C
  38 C Body
  39 C
  40 C Set LPRINT=.TRUE. for debugging
  41       dwa16=2.0d0**(1.0d0/6.0d0)
  42       lprint=.false.
  43       itypro=20
  44 C Assign virtual-bond length
  45       vbl=3.8D0
  46       vblinv=1.0D0/vbl
  47       vblinv2=vblinv*vblinv
  48       call card_concat(controlcard,.true.)
  49       wname(4)="WCORRH"
  50       do i=1,n_ene
  51         key = wname(i)(:ilen(wname(i)))
  52         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
  53       enddo
  54
  55       write (iout,*) "iparm",iparm," myparm",myparm
  56 c If reading not own parameters, skip assignment
  57
  58       if (iparm.eq.myparm .or. .not.separate_parset) then
  59
  60 c
  61 c Setup weights for UNRES
  62 c
  63       wsc=ww(1)
  64       wscp=ww(2)
  65       welec=ww(3)
  66       wcorr=ww(4)
  67       wcorr5=ww(5)
  68       wcorr6=ww(6)
  69       wel_loc=ww(7)
  70       wturn3=ww(8)
  71       wturn4=ww(9)
  72       wturn6=ww(10)
  73       wang=ww(11)
  74       wscloc=ww(12)
  75       wtor=ww(13)
  76       wtor_d=ww(14)
  77       wvdwpp=ww(16)
  78       wstrain=ww(15)
  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.onelett(i) .or. t2.ne.onelett(j)
 600      &        .or. t3.ne.onelett(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       do i=-ntyp,-1
 687         isccortyp(i)=-isccortyp(-i)
 688        enddo
 689        iscprol=isccortyp(20)
 690 c      write (iout,*) 'ntortyp',ntortyp
 691       maxinter=3
 692 cc maxinter is maximum interaction sites
 693       do l=1,maxinter
 694       do i=1,nsccortyp
 695         do j=1,nsccortyp
 696           read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
 697           v0ijsccor=0.0d0
 698           v0ijsccor1=0.0d0
 699           v0ijsccor2=0.0d0
 700           v0ijsccor3=0.0d0
 701           si=-1.0d0
 702           nterm_sccor(-i,j)=nterm_sccor(i,j)
 703           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 704           nterm_sccor(i,-j)=nterm_sccor(i,j)
 705           do k=1,nterm_sccor(i,j)
 706             read (isccor,*) kk,v1sccor(k,l,i,j)
 707      &    ,v2sccor(k,l,i,j)
 708              if (j.eq.iscprol) then
 709              if (i.eq.isccortyp(10)) then
 710              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 711              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 712              else
 713               v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 714      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 715              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 716      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 717              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 718              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 719              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 720              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 721              endif
 722             else
 723              if (i.eq.isccortyp(10)) then
 724              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 725              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 726              else
 727                if (j.eq.isccortyp(10)) then
 728              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 729              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 730                else
 731              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 732              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 733              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 734              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 735              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 736              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 737                 endif
 738                endif
 739              endif
 740             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 741             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 742             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 743             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 744             si=-si
 745           enddo
 746           do k=1,nlor_sccor(i,j)
 747             read (isccor,*) kk,vlor1sccor(k,i,j),
 748      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 749             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 750      &(1+vlor3sccor(k,i,j)**2)
 751           enddo
 752           v0sccor(l,i,j)=v0ijsccor
 753           v0sccor(l,-i,j)=v0ijsccor1
 754           v0sccor(l,i,-j)=v0ijsccor2
 755           v0sccor(l,-i,-j)=v0ijsccor3
 756         enddo
 757       enddo
 758       enddo
 759       close (isccor)
 760
 761       if (lprint) then
 762         write (iout,'(/a/)') 'Torsional constants:'
 763         do i=1,nsccortyp
 764           do j=1,nsccortyp
 765             write (iout,*) 'ityp',i,' jtyp',j
 766             write (iout,*) 'Fourier constants'
 767             do k=1,nterm_sccor(i,j)
 768            write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 769             enddo
 770             write (iout,*) 'Lorenz constants'
 771             do k=1,nlor_sccor(i,j)
 772               write (iout,'(3(1pe15.5))')
 773      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 774             enddo
 775           enddo
 776         enddo
 777       endif
 778
 779 C
 780 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 781 C         interaction energy of the Gly, Ala, and Pro prototypes.
 782 C
 783       read (ifourier,*) nloctyp
 784       do i=0,nloctyp-1
 785         read (ifourier,*)
 786         read (ifourier,*) (b(ii),ii=1,13)
 787         if (lprint) then
 788         write (iout,*) 'Type',i
 789         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 790         endif
 791         B1(1,i)  = b(3)
 792         B1(2,i)  = b(5)
 793         B1(1,-i) = b(3)
 794         B1(2,-i) = -b(5)
 795 c        b1(1,i)=0.0d0
 796 c        b1(2,i)=0.0d0
 797         B1tilde(1,i) = b(3)
 798         B1tilde(2,i) =-b(5)
 799         B1tilde(1,-i) =-b(3)
 800         B1tilde(2,-i) =b(5)
 801 c        b1tilde(1,i)=0.0d0
 802 c        b1tilde(2,i)=0.0d0
 803         B2(1,i)  = b(2)
 804         B2(2,i)  = b(4)
 805         B2(1,-i)  =b(2)
 806         B2(2,-i)  =-b(4)
 807
 808 c        b2(1,i)=0.0d0
 809 c        b2(2,i)=0.0d0
 810         CC(1,1,i)= b(7)
 811         CC(2,2,i)=-b(7)
 812         CC(2,1,i)= b(9)
 813         CC(1,2,i)= b(9)
 814         CC(1,1,-i)= b(7)
 815         CC(2,2,-i)=-b(7)
 816         CC(2,1,-i)=-b(9)
 817         CC(1,2,-i)=-b(9)
 818 c        CC(1,1,i)=0.0d0
 819 c        CC(2,2,i)=0.0d0
 820 c        CC(2,1,i)=0.0d0
 821 c        CC(1,2,i)=0.0d0
 822         Ctilde(1,1,i)=b(7)
 823         Ctilde(1,2,i)=b(9)
 824         Ctilde(2,1,i)=-b(9)
 825         Ctilde(2,2,i)=b(7)
 826         Ctilde(1,1,-i)=b(7)
 827         Ctilde(1,2,-i)=-b(9)
 828         Ctilde(2,1,-i)=b(9)
 829         Ctilde(2,2,-i)=b(7)
 830
 831 c        Ctilde(1,1,i)=0.0d0
 832 c        Ctilde(1,2,i)=0.0d0
 833 c        Ctilde(2,1,i)=0.0d0
 834 c        Ctilde(2,2,i)=0.0d0
 835         DD(1,1,i)= b(6)
 836         DD(2,2,i)=-b(6)
 837         DD(2,1,i)= b(8)
 838         DD(1,2,i)= b(8)
 839         DD(1,1,-i)= b(6)
 840         DD(2,2,-i)=-b(6)
 841         DD(2,1,-i)=-b(8)
 842         DD(1,2,-i)=-b(8)
 843 c        DD(1,1,i)=0.0d0
 844 c        DD(2,2,i)=0.0d0
 845 c        DD(2,1,i)=0.0d0
 846 c        DD(1,2,i)=0.0d0
 847         Dtilde(1,1,i)=b(6)
 848         Dtilde(1,2,i)=b(8)
 849         Dtilde(2,1,i)=-b(8)
 850         Dtilde(2,2,i)=b(6)
 851         Dtilde(1,1,-i)=b(6)
 852         Dtilde(1,2,-i)=-b(8)
 853         Dtilde(2,1,-i)=b(8)
 854         Dtilde(2,2,-i)=b(6)
 855
 856 c        Dtilde(1,1,i)=0.0d0
 857 c        Dtilde(1,2,i)=0.0d0
 858 c        Dtilde(2,1,i)=0.0d0
 859 c        Dtilde(2,2,i)=0.0d0
 860         EE(1,1,i)= b(10)+b(11)
 861         EE(2,2,i)=-b(10)+b(11)
 862         EE(2,1,i)= b(12)-b(13)
 863         EE(1,2,i)= b(12)+b(13)
 864         EE(1,1,-i)= b(10)+b(11)
 865         EE(2,2,-i)=-b(10)+b(11)
 866         EE(2,1,-i)=-b(12)+b(13)
 867         EE(1,2,-i)=-b(12)-b(13)
 868
 869 c        ee(1,1,i)=1.0d0
 870 c        ee(2,2,i)=1.0d0
 871 c        ee(2,1,i)=0.0d0
 872 c        ee(1,2,i)=0.0d0
 873 c        ee(2,1,i)=ee(1,2,i)
 874       enddo
 875       if (lprint) then
 876       do i=1,nloctyp
 877         write (iout,*) 'Type',i
 878         write (iout,*) 'B1'
 879 c        write (iout,'(f10.5)') B1(:,i)
 880         write(iout,*) B1(1,i),B1(2,i)
 881         write (iout,*) 'B2'
 882 c        write (iout,'(f10.5)') B2(:,i)
 883         write(iout,*) B2(1,i),B2(2,i)
 884         write (iout,*) 'CC'
 885         do j=1,2
 886           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 887         enddo
 888         write(iout,*) 'DD'
 889         do j=1,2
 890           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 891         enddo
 892         write(iout,*) 'EE'
 893         do j=1,2
 894           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 895         enddo
 896       enddo
 897       endif
 898 C
 899 C Read electrostatic-interaction parameters
 900 C
 901       if (lprint) then
 902         write (iout,'(/a)') 'Electrostatic interaction constants:'
 903         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 904      &            'IT','JT','APP','BPP','AEL6','AEL3'
 905       endif
 906       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 907       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 908       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 909       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 910       close (ielep)
 911       do i=1,2
 912         do j=1,2
 913         rri=rpp(i,j)**6
 914         app (i,j)=epp(i,j)*rri*rri
 915         bpp (i,j)=-2.0D0*epp(i,j)*rri
 916         ael6(i,j)=elpp6(i,j)*4.2D0**6
 917         ael3(i,j)=elpp3(i,j)*4.2D0**3
 918         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 919      &                    ael6(i,j),ael3(i,j)
 920         enddo
 921       enddo
 922 C
 923 C Read side-chain interaction parameters.
 924 C
 925       read (isidep,*) ipot,expon
 926       if (ipot.lt.1 .or. ipot.gt.5) then
 927         write (iout,'(2a)') 'Error while reading SC interaction',
 928      &               'potential file - unknown potential type.'
 929         stop
 930       endif
 931       expon2=expon/2
 932       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 933      & ', exponents are ',expon,2*expon
 934       goto (10,20,30,30,40) ipot
 935 C----------------------- LJ potential ---------------------------------
 936    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 937       if (lprint) then
 938         write (iout,'(/a/)') 'Parameters of the LJ 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,a)') 'residue','sigma'
 943         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 944       endif
 945       goto 50
 946 C----------------------- LJK potential --------------------------------
 947    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 948      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 949       if (lprint) then
 950         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 951         write (iout,'(a/)') 'The epsilon array:'
 952         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 953         write (iout,'(/a)') 'One-body parameters:'
 954         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 955         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 956      &        i=1,ntyp)
 957       endif
 958       goto 50
 959 C---------------------- GB or BP potential -----------------------------
 960    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 961      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 962      &  (alp(i),i=1,ntyp)
 963 C For the GB potential convert sigma'**2 into chi'
 964       if (ipot.eq.4) then
 965         do i=1,ntyp
 966           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 967         enddo
 968       endif
 969       if (lprint) then
 970         write (iout,'(/a/)') 'Parameters of the BP potential:'
 971         write (iout,'(a/)') 'The epsilon array:'
 972         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 973         write (iout,'(/a)') 'One-body parameters:'
 974         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 975      &       '    chip  ','    alph  '
 976         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 977      &                     chip(i),alp(i),i=1,ntyp)
 978       endif
 979       goto 50
 980 C--------------------- GBV potential -----------------------------------
 981    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 982      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 983      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 984       if (lprint) then
 985         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 986         write (iout,'(a/)') 'The epsilon array:'
 987         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 988         write (iout,'(/a)') 'One-body parameters:'
 989         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 990      &      's||/s_|_^2','    chip  ','    alph  '
 991         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 992      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 993       endif
 994    50 continue
 995       close (isidep)
 996 C-----------------------------------------------------------------------
 997 C Calculate the "working" parameters of SC interactions.
 998       do i=2,ntyp
 999         do j=1,i-1
1000           eps(i,j)=eps(j,i)
1001         enddo
1002       enddo
1003       do i=1,ntyp
1004         do j=i,ntyp
1005           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1006           sigma(j,i)=sigma(i,j)
1007           rs0(i,j)=dwa16*sigma(i,j)
1008           rs0(j,i)=rs0(i,j)
1009         enddo
1010       enddo
1011       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1012      & 'Working parameters of the SC interactions:',
1013      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1014      & '  chi1   ','   chi2   '
1015       do i=1,ntyp
1016         do j=i,ntyp
1017           epsij=eps(i,j)
1018           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1019             rrij=sigma(i,j)
1020           else
1021             rrij=rr0(i)+rr0(j)
1022           endif
1023           r0(i,j)=rrij
1024           r0(j,i)=rrij
1025           rrij=rrij**expon
1026           epsij=eps(i,j)
1027           sigeps=dsign(1.0D0,epsij)
1028           epsij=dabs(epsij)
1029           aa(i,j)=epsij*rrij*rrij
1030           bb(i,j)=-sigeps*epsij*rrij
1031           aa(j,i)=aa(i,j)
1032           bb(j,i)=bb(i,j)
1033           if (ipot.gt.2) then
1034             sigt1sq=sigma0(i)**2
1035             sigt2sq=sigma0(j)**2
1036             sigii1=sigii(i)
1037             sigii2=sigii(j)
1038             ratsig1=sigt2sq/sigt1sq
1039             ratsig2=1.0D0/ratsig1
1040             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1041             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1042             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1043           else
1044             rsum_max=sigma(i,j)
1045           endif
1046 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1047             sigmaii(i,j)=rsum_max
1048             sigmaii(j,i)=rsum_max
1049 c         else
1050 c           sigmaii(i,j)=r0(i,j)
1051 c           sigmaii(j,i)=r0(i,j)
1052 c         endif
1053 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1054           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1055             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1056             augm(i,j)=epsij*r_augm**(2*expon)
1057 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1058             augm(j,i)=augm(i,j)
1059           else
1060             augm(i,j)=0.0D0
1061             augm(j,i)=0.0D0
1062           endif
1063           if (lprint) then
1064             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1065      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1066      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1067           endif
1068         enddo
1069       enddo
1070 C
1071 C Define the SC-p interaction constants
1072 C
1073 #ifdef OLDSCP
1074       do i=1,20
1075 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1076 C helix formation)
1077 c       aad(i,1)=0.3D0*4.0D0**12
1078 C Following line for constants currently implemented
1079 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1080         aad(i,1)=1.5D0*4.0D0**12
1081 c       aad(i,1)=0.17D0*5.6D0**12
1082         aad(i,2)=aad(i,1)
1083 C "Soft" SC-p repulsion
1084         bad(i,1)=0.0D0
1085 C Following line for constants currently implemented
1086 c       aad(i,1)=0.3D0*4.0D0**6
1087 C "Hard" SC-p repulsion
1088         bad(i,1)=3.0D0*4.0D0**6
1089 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1090         bad(i,2)=bad(i,1)
1091 c       aad(i,1)=0.0D0
1092 c       aad(i,2)=0.0D0
1093 c       bad(i,1)=1228.8D0
1094 c       bad(i,2)=1228.8D0
1095       enddo
1096 #else
1097 C
1098 C 8/9/01 Read the SC-p interaction constants from file
1099 C
1100       do i=1,ntyp
1101         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1102       enddo
1103       do i=1,ntyp
1104         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1105         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1106         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1107         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1108       enddo
1109
1110       if (lprint) then
1111         write (iout,*) "Parameters of SC-p interactions:"
1112         do i=1,20
1113           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1114      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1115         enddo
1116       endif
1117 #endif
1118 C
1119 C Define the constants of the disulfide bridge
1120 C
1121       ebr=-5.50D0
1122 c
1123 c Old arbitrary potential - commented out.
1124 c
1125 c      dbr= 4.20D0
1126 c      fbr= 3.30D0
1127 c
1128 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1129 c energy surface of diethyl disulfide.
1130 c A. Liwo and U. Kozlowska, 11/24/03
1131 c
1132       D0CM = 3.78d0
1133       AKCM = 15.1d0
1134       AKTH = 11.0d0
1135       AKCT = 12.0d0
1136       V1SS =-1.08d0
1137       V2SS = 7.61d0
1138       V3SS = 13.7d0
1139
1140       if (lprint) then
1141       write (iout,'(/a)') "Disulfide bridge parameters:"
1142       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1143       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1144       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1145       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1146      & ' v3ss:',v3ss
1147       endif
1148       return
1149       end