source/cluster/wham/src-M/parmread.F

   1       subroutine parmread
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8       include 'sizesclu.dat'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.INTERACT'
  12       include 'COMMON.GEO'
  13       include 'COMMON.LOCAL'
  14       include 'COMMON.TORSION'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.NAMES'
  17       include 'COMMON.SBRIDGE'
  18       include 'COMMON.SCCOR'
  19       include 'COMMON.SCROT'
  20       character*1 t1,t2,t3
  21       character*1 onelett(4) /"G","A","P","D"/
  22       logical lprint
  23       dimension blower(3,3,maxlob)
  24       double precision ip,mp
  25 C
  26 C Body
  27 C
  28 C Set LPRINT=.TRUE. for debugging
  29       dwa16=2.0d0**(1.0d0/6.0d0)
  30       lprint=.false.
  31       itypro=20
  32 C Assign virtual-bond length
  33       vbl=3.8D0
  34       vblinv=1.0D0/vbl
  35       vblinv2=vblinv*vblinv
  36 #ifdef CRYST_BOND
  37       read (ibond,*) vbldp0,akp
  38       do i=1,ntyp
  39         nbondterm(i)=1
  40         read (ibond,*) vbldsc0(1,i),aksc(1,i)
  41         dsc(i) = vbldsc0(1,i)
  42         if (i.eq.10) then
  43           dsc_inv(i)=0.0D0
  44         else
  45           dsc_inv(i)=1.0D0/dsc(i)
  46         endif
  47       enddo
  48 #else
  49       read (ibond,*) ijunk,vbldp0,akp,rjunk
  50       do i=1,ntyp
  51         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
  52      &   j=1,nbondterm(i))
  53         dsc(i) = vbldsc0(1,i)
  54         if (i.eq.10) then
  55           dsc_inv(i)=0.0D0
  56         else
  57           dsc_inv(i)=1.0D0/dsc(i)
  58         endif
  59       enddo
  60 #endif
  61       if (lprint) then
  62         write(iout,'(/a/)')"Force constants virtual bonds:"
  63         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
  64      &   'inertia','Pstok'
  65         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
  66         do i=1,ntyp
  67           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
  68      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
  69           do j=2,nbondterm(i)
  70             write (iout,'(13x,3f10.5)')
  71      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
  72           enddo
  73         enddo
  74       endif
  75 #ifdef CRYST_THETA
  76 C
  77 C Read the parameters of the probability distribution/energy expression
  78 C of the virtual-bond valence angles theta
  79 C
  80       do i=1,ntyp
  81         read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
  82         read (ithep,*) (polthet(j,i),j=0,3)
  83         read (ithep,*) (gthet(j,i),j=1,3)
  84         read (ithep,*) theta0(i),sig0(i),sigc0(i)
  85         sigc0(i)=sigc0(i)**2
  86       enddo
  87       close (ithep)
  88       if (lprint) then
  89 c       write (iout,'(a)')
  90 c    &    'Parameters of the virtual-bond valence angles:'
  91 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
  92 c    & '    ATHETA0   ','         A1   ','        A2    ',
  93 c    & '        B1    ','         B2   '
  94 c       do i=1,ntyp
  95 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
  96 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
  97 c       enddo
  98 c       write (iout,'(/a/9x,5a/79(1h-))')
  99 c    & 'Parameters of the expression for sigma(theta_c):',
 100 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 101 c    & '     ALPH3    ','    SIGMA0C   '
 102 c       do i=1,ntyp
 103 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 104 c    &      (polthet(j,i),j=0,3),sigc0(i)
 105 c       enddo
 106 c       write (iout,'(/a/9x,5a/79(1h-))')
 107 c    & 'Parameters of the second gaussian:',
 108 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 109 c    & '        G2    ','         G3   '
 110 c       do i=1,ntyp
 111 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 112 c    &       sig0(i),(gthet(j,i),j=1,3)
 113 c       enddo
 114         write (iout,'(a)')
 115      &    'Parameters of the virtual-bond valence angles:'
 116         write (iout,'(/a/9x,5a/79(1h-))')
 117      & 'Coefficients of expansion',
 118      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 119      & '   b1*10^1    ','    b2*10^1   '
 120         do i=1,ntyp
 121           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 122      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
 123         enddo
 124         write (iout,'(/a/9x,5a/79(1h-))')
 125      & 'Parameters of the expression for sigma(theta_c):',
 126      & ' alpha0       ','  alph1       ',' alph2        ',
 127      & ' alhp3        ','   sigma0c    '
 128         do i=1,ntyp
 129           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 130      &      (polthet(j,i),j=0,3),sigc0(i)
 131         enddo
 132         write (iout,'(/a/9x,5a/79(1h-))')
 133      & 'Parameters of the second gaussian:',
 134      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 135      & '        G2    ','   G3*10^1    '
 136         do i=1,ntyp
 137           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 138      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 139         enddo
 140       endif
 141 #else
 142 C
 143 C Read the parameters of Utheta determined from ab initio surfaces
 144 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 145 C
 146       read (ithep,*) nthetyp,ntheterm,ntheterm2,
 147      &  ntheterm3,nsingle,ndouble
 148       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 149       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 150       do i=1,maxthetyp
 151         do j=1,maxthetyp
 152           do k=1,maxthetyp
 153             aa0thet(i,j,k)=0.0d0
 154             do l=1,ntheterm
 155               aathet(l,i,j,k)=0.0d0
 156             enddo
 157             do l=1,ntheterm2
 158               do m=1,nsingle
 159                 bbthet(m,l,i,j,k)=0.0d0
 160                 ccthet(m,l,i,j,k)=0.0d0
 161                 ddthet(m,l,i,j,k)=0.0d0
 162                 eethet(m,l,i,j,k)=0.0d0
 163               enddo
 164             enddo
 165             do l=1,ntheterm3
 166               do m=1,ndouble
 167                 do mm=1,ndouble
 168                  ffthet(mm,m,l,i,j,k)=0.0d0
 169                  ggthet(mm,m,l,i,j,k)=0.0d0
 170                 enddo
 171               enddo
 172             enddo
 173           enddo
 174         enddo
 175       enddo
 176       do i=1,nthetyp
 177         do j=1,nthetyp
 178           do k=1,nthetyp
 179             read (ithep,'(3a)') res1,res2,res3
 180             read (ithep,*) aa0thet(i,j,k)
 181             read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
 182             read (ithep,*)
 183      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
 184      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
 185      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
 186      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
 187             read (ithep,*)
 188      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
 189      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
 190      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 191           enddo
 192         enddo
 193       enddo
 194 C
 195 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 196 C coefficients of theta-and-gamma-dependent terms are zero.
 197 C
 198       do i=1,nthetyp
 199         do j=1,nthetyp
 200           do l=1,ntheterm
 201             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
 202             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
 203           enddo
 204           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
 205           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
 206         enddo
 207         do l=1,ntheterm
 208           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
 209         enddo
 210         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
 211       enddo
 212 C
 213 C Control printout of the coefficients of virtual-bond-angle potentials
 214 C
 215       if (lprint) then
 216         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 217         do i=1,nthetyp+1
 218           do j=1,nthetyp+1
 219             do k=1,nthetyp+1
 220               write (iout,'(//4a)')
 221      &         'Type ',onelett(i),onelett(j),onelett(k)
 222               write (iout,'(//a,10x,a)') " l","a[l]"
 223               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
 224               write (iout,'(i2,1pe15.5)')
 225      &           (l,aathet(l,i,j,k),l=1,ntheterm)
 226             do l=1,ntheterm2
 227               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 228      &          "b",l,"c",l,"d",l,"e",l
 229               do m=1,nsingle
 230                 write (iout,'(i2,4(1pe15.5))') m,
 231      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
 232      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
 233               enddo
 234             enddo
 235             do l=1,ntheterm3
 236               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 237      &          "f+",l,"f-",l,"g+",l,"g-",l
 238               do m=2,ndouble
 239                 do n=1,m-1
 240                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 241      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
 242      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
 243                 enddo
 244               enddo
 245             enddo
 246           enddo
 247         enddo
 248       enddo
 249       call flush(iout)
 250       endif
 251 #endif
 252
 253 #ifdef CRYST_SC
 254 C
 255 C Read the parameters of the probability distribution/energy expression
 256 C of the side chains.
 257 C
 258       do i=1,ntyp
 259         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 260         if (i.eq.10) then
 261           dsc_inv(i)=0.0D0
 262         else
 263           dsc_inv(i)=1.0D0/dsc(i)
 264         endif
 265         if (i.ne.10) then
 266         do j=1,nlob(i)
 267           do k=1,3
 268             do l=1,3
 269               blower(l,k,j)=0.0D0
 270             enddo
 271           enddo
 272         enddo
 273         bsc(1,i)=0.0D0
 274         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 275         do j=2,nlob(i)
 276           read (irotam,*) bsc(j,i)
 277           read (irotam,*) (censc(k,j,i),k=1,3),
 278      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 279         enddo
 280         do j=1,nlob(i)
 281           do k=1,3
 282             do l=1,k
 283               akl=0.0D0
 284               do m=1,3
 285                 akl=akl+blower(k,m,j)*blower(l,m,j)
 286               enddo
 287               gaussc(k,l,j,i)=akl
 288               gaussc(l,k,j,i)=akl
 289             enddo
 290           enddo
 291         enddo
 292         endif
 293       enddo
 294       close (irotam)
 295       if (lprint) then
 296         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 297         do i=1,ntyp
 298           nlobi=nlob(i)
 299           if (nlobi.gt.0) then
 300           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 301      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 302 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 303 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 304 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 305 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 306            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 307      &                             'log h',(bsc(j,i),j=1,nlobi)
 308            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 309      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 310 c          write (iout,'(a)')
 311 c         do j=1,nlobi
 312 c           ind=0
 313 c           do k=1,3
 314 c             do l=1,k
 315 c              ind=ind+1
 316 c              blower(k,l,j)=gaussc(ind,j,i)
 317 c             enddo
 318 c           enddo
 319 c         enddo
 320           do k=1,3
 321             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 322      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 323           enddo
 324           endif
 325         enddo
 326       endif
 327 #else
 328 C
 329 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 330 C added by Urszula Kozlowska 07/11/2007
 331 C
 332       do i=1,ntyp
 333         read (irotam,*)
 334        if (i.eq.10) then
 335          read (irotam,*)
 336        else
 337          do j=1,65
 338            read(irotam,*) sc_parmin(j,i)
 339          enddo
 340        endif
 341       enddo
 342 #endif
 343       close(irotam)
 344 #ifdef CRYST_TOR
 345 C
 346 C Read torsional parameters in old format
 347 C
 348       read (itorp,*) ntortyp,nterm_old
 349       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 350       read (itorp,*) (itortyp(i),i=1,ntyp)
 351       do i=1,ntortyp
 352         do j=1,ntortyp
 353           read (itorp,'(a)')
 354           do k=1,nterm_old
 355             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 356           enddo
 357         enddo
 358       enddo
 359       close (itorp)
 360       if (lprint) then
 361         write (iout,'(/a/)') 'Torsional constants:'
 362         do i=1,ntortyp
 363           do j=1,ntortyp
 364             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 365             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 366           enddo
 367         enddo
 368       endif
 369 #else
 370 C
 371 C Read torsional parameters
 372 C
 373       read (itorp,*) ntortyp
 374       read (itorp,*) (itortyp(i),i=1,ntyp)
 375       write (iout,*) 'ntortyp',ntortyp
 376       do i=1,ntortyp
 377         do j=1,ntortyp
 378           read (itorp,*) nterm(i,j),nlor(i,j)
 379           v0ij=0.0d0
 380           si=-1.0d0
 381           do k=1,nterm(i,j)
 382             read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
 383             v0ij=v0ij+si*v1(k,i,j)
 384             si=-si
 385           enddo
 386           do k=1,nlor(i,j)
 387             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 388             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 389           enddo
 390           v0(i,j)=v0ij
 391         enddo
 392       enddo
 393       close (itorp)
 394       if (lprint) then
 395         write (iout,'(/a/)') 'Torsional constants:'
 396         do i=1,ntortyp
 397           do j=1,ntortyp
 398             write (iout,*) 'ityp',i,' jtyp',j
 399             write (iout,*) 'Fourier constants'
 400             do k=1,nterm(i,j)
 401               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 402             enddo
 403             write (iout,*) 'Lorenz constants'
 404             do k=1,nlor(i,j)
 405               write (iout,'(3(1pe15.5))')
 406      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 407             enddo
 408           enddo
 409         enddo
 410       endif
 411 C
 412 C 6/23/01 Read parameters for double torsionals
 413 C
 414       do i=1,ntortyp
 415         do j=1,ntortyp
 416           do k=1,ntortyp
 417             read (itordp,'(3a1)') t1,t2,t3
 418             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 419      &        .or. t3.ne.onelett(k)) then
 420               write (iout,*) "Error in double torsional parameter file",
 421      &         i,j,k,t1,t2,t3
 422                stop "Error in double torsional parameter file"
 423             endif
 424             read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
 425             read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
 426             read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
 427             read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
 428             read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
 429             read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
 430      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
 431           enddo
 432         enddo
 433       enddo
 434       if (lprint) then
 435       write (iout,*)
 436       write (iout,*) 'Constants for double torsionals'
 437       do i=1,ntortyp
 438         do j=1,ntortyp
 439           do k=1,ntortyp
 440             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 441      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 442             write (iout,*)
 443             write (iout,*) 'Single angles:'
 444             do l=1,ntermd_1(i,j,k)
 445               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 446      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 447      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 448             enddo
 449             write (iout,*)
 450             write (iout,*) 'Pairs of angles:'
 451             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 452             do l=1,ntermd_2(i,j,k)
 453               write (iout,'(i5,20f10.5)')
 454      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 455             enddo
 456             write (iout,*)
 457             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 458             do l=1,ntermd_2(i,j,k)
 459               write (iout,'(i5,20f10.5)')
 460      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 461             enddo
 462             write (iout,*)
 463           enddo
 464         enddo
 465       enddo
 466       endif
 467 #endif
 468 C
 469 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 470 C         correlation energies.
 471 C
 472       read (isccor,*) nterm_sccor
 473       do i=1,20
 474         do j=1,20
 475           read (isccor,'(a)')
 476           do k=1,nterm_sccor
 477             read (isccor,*)
 478      &        kk,v1sccor(k,i,j),v2sccor(k,i,j)
 479           enddo
 480         enddo
 481       enddo
 482       close (isccor)
 483       if (lprint) then
 484         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 485         do i=1,20
 486           do j=1,20
 487             write (iout,*) 'ityp',i,' jtyp',j
 488             do k=1,nterm_sccor
 489               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
 490             enddo
 491           enddo
 492         enddo
 493       endif
 494 C
 495 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 496 C         interaction energy of the Gly, Ala, and Pro prototypes.
 497 C
 498       read (ifourier,*) nloctyp
 499       do i=1,nloctyp
 500         read (ifourier,*)
 501         read (ifourier,*) (b(ii,i),ii=1,13)
 502         if (lprint) then
 503         write (iout,*) 'Type',i
 504         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 505         endif
 506         B1(1,i)  = b(3,i)
 507         B1(2,i)  = b(5,i)
 508         B1tilde(1,i) = b(3,i)
 509         B1tilde(2,i) =-b(5,i)
 510         B2(1,i)  = b(2,i)
 511         B2(2,i)  = b(4,i)
 512         CC(1,1,i)= b(7,i)
 513         CC(2,2,i)=-b(7,i)
 514         CC(2,1,i)= b(9,i)
 515         CC(1,2,i)= b(9,i)
 516         Ctilde(1,1,i)=b(7,i)
 517         Ctilde(1,2,i)=b(9,i)
 518         Ctilde(2,1,i)=-b(9,i)
 519         Ctilde(2,2,i)=b(7,i)
 520         DD(1,1,i)= b(6,i)
 521         DD(2,2,i)=-b(6,i)
 522         DD(2,1,i)= b(8,i)
 523         DD(1,2,i)= b(8,i)
 524         Dtilde(1,1,i)=b(6,i)
 525         Dtilde(1,2,i)=b(8,i)
 526         Dtilde(2,1,i)=-b(8,i)
 527         Dtilde(2,2,i)=b(6,i)
 528         EE(1,1,i)= b(10,i)+b(11,i)
 529         EE(2,2,i)=-b(10,i)+b(11,i)
 530         EE(2,1,i)= b(12,i)-b(13,i)
 531         EE(1,2,i)= b(12,i)+b(13,i)
 532       enddo
 533       if (lprint) then
 534       do i=1,nloctyp
 535         write (iout,*) 'Type',i
 536         write (iout,*) 'B1'
 537 c        write (iout,'(f10.5)') B1(:,i)
 538         write(iout,*) B1(1,i),B1(2,i)
 539         write (iout,*) 'B2'
 540 c        write (iout,'(f10.5)') B2(:,i)
 541         write(iout,*) B2(1,i),B2(2,i)
 542         write (iout,*) 'CC'
 543         do j=1,2
 544           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 545         enddo
 546         write(iout,*) 'DD'
 547         do j=1,2
 548           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 549         enddo
 550         write(iout,*) 'EE'
 551         do j=1,2
 552           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 553         enddo
 554       enddo
 555       endif
 556 C
 557 C Read electrostatic-interaction parameters
 558 C
 559       if (lprint) then
 560         write (iout,'(/a)') 'Electrostatic interaction constants:'
 561         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 562      &            'IT','JT','APP','BPP','AEL6','AEL3'
 563       endif
 564       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 565       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 566       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 567       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 568       close (ielep)
 569       do i=1,2
 570         do j=1,2
 571         rri=rpp(i,j)**6
 572         app (i,j)=epp(i,j)*rri*rri
 573         bpp (i,j)=-2.0D0*epp(i,j)*rri
 574         ael6(i,j)=elpp6(i,j)*4.2D0**6
 575         ael3(i,j)=elpp3(i,j)*4.2D0**3
 576         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 577      &                    ael6(i,j),ael3(i,j)
 578         enddo
 579       enddo
 580 C
 581 C Read side-chain interaction parameters.
 582 C
 583       read (isidep,*) ipot,expon
 584       if (ipot.lt.1 .or. ipot.gt.5) then
 585         write (iout,'(2a)') 'Error while reading SC interaction',
 586      &               'potential file - unknown potential type.'
 587         stop
 588       endif
 589       expon2=expon/2
 590       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 591      & ', exponents are ',expon,2*expon
 592       goto (10,20,30,30,40) ipot
 593 C----------------------- LJ potential ---------------------------------
 594    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 595       if (lprint) then
 596         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 597         write (iout,'(a/)') 'The epsilon array:'
 598         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 599         write (iout,'(/a)') 'One-body parameters:'
 600         write (iout,'(a,4x,a)') 'residue','sigma'
 601         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 602       endif
 603       goto 50
 604 C----------------------- LJK potential --------------------------------
 605    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 606      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 607       if (lprint) then
 608         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 609         write (iout,'(a/)') 'The epsilon array:'
 610         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 611         write (iout,'(/a)') 'One-body parameters:'
 612         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 613         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 614      &        i=1,ntyp)
 615       endif
 616       goto 50
 617 C---------------------- GB or BP potential -----------------------------
 618    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 619      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 620      &  (alp(i),i=1,ntyp)
 621 C For the GB potential convert sigma'**2 into chi'
 622       if (ipot.eq.4) then
 623         do i=1,ntyp
 624           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 625         enddo
 626       endif
 627       if (lprint) then
 628         write (iout,'(/a/)') 'Parameters of the BP potential:'
 629         write (iout,'(a/)') 'The epsilon array:'
 630         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 631         write (iout,'(/a)') 'One-body parameters:'
 632         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 633      &       '    chip  ','    alph  '
 634         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 635      &                     chip(i),alp(i),i=1,ntyp)
 636       endif
 637       goto 50
 638 C--------------------- GBV potential -----------------------------------
 639    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 640      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 641      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 642       if (lprint) then
 643         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 644         write (iout,'(a/)') 'The epsilon array:'
 645         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 646         write (iout,'(/a)') 'One-body parameters:'
 647         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 648      &      's||/s_|_^2','    chip  ','    alph  '
 649         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 650      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 651       endif
 652    50 continue
 653       close (isidep)
 654 C-----------------------------------------------------------------------
 655 C Calculate the "working" parameters of SC interactions.
 656       do i=2,ntyp
 657         do j=1,i-1
 658           eps(i,j)=eps(j,i)
 659         enddo
 660       enddo
 661       do i=1,ntyp
 662         do j=i,ntyp
 663           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 664           sigma(j,i)=sigma(i,j)
 665           rs0(i,j)=dwa16*sigma(i,j)
 666           rs0(j,i)=rs0(i,j)
 667         enddo
 668       enddo
 669       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 670      & 'Working parameters of the SC interactions:',
 671      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 672      & '  chi1   ','   chi2   '
 673       do i=1,ntyp
 674         do j=i,ntyp
 675           epsij=eps(i,j)
 676           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 677             rrij=sigma(i,j)
 678           else
 679             rrij=rr0(i)+rr0(j)
 680           endif
 681           r0(i,j)=rrij
 682           r0(j,i)=rrij
 683           rrij=rrij**expon
 684           epsij=eps(i,j)
 685           sigeps=dsign(1.0D0,epsij)
 686           epsij=dabs(epsij)
 687           aa(i,j)=epsij*rrij*rrij
 688           bb(i,j)=-sigeps*epsij*rrij
 689           aa(j,i)=aa(i,j)
 690           bb(j,i)=bb(i,j)
 691           if (ipot.gt.2) then
 692             sigt1sq=sigma0(i)**2
 693             sigt2sq=sigma0(j)**2
 694             sigii1=sigii(i)
 695             sigii2=sigii(j)
 696             ratsig1=sigt2sq/sigt1sq
 697             ratsig2=1.0D0/ratsig1
 698             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 699             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 700             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 701           else
 702             rsum_max=sigma(i,j)
 703           endif
 704 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 705             sigmaii(i,j)=rsum_max
 706             sigmaii(j,i)=rsum_max
 707 c         else
 708 c           sigmaii(i,j)=r0(i,j)
 709 c           sigmaii(j,i)=r0(i,j)
 710 c         endif
 711 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 712           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 713             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 714             augm(i,j)=epsij*r_augm**(2*expon)
 715 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 716             augm(j,i)=augm(i,j)
 717           else
 718             augm(i,j)=0.0D0
 719             augm(j,i)=0.0D0
 720           endif
 721           if (lprint) then
 722             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 723      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 724      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 725           endif
 726         enddo
 727       enddo
 728 C
 729 C Define the SC-p interaction constants
 730 C
 731       do i=1,20
 732         do j=1,2
 733           eps_scp(i,j)=-1.5d0
 734           rscp(i,j)=4.0d0
 735         enddo
 736       enddo
 737 #ifdef OLDSCP
 738       do i=1,20
 739 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 740 C helix formation)
 741 c       aad(i,1)=0.3D0*4.0D0**12
 742 C Following line for constants currently implemented
 743 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 744         aad(i,1)=1.5D0*4.0D0**12
 745 c       aad(i,1)=0.17D0*5.6D0**12
 746         aad(i,2)=aad(i,1)
 747 C "Soft" SC-p repulsion
 748         bad(i,1)=0.0D0
 749 C Following line for constants currently implemented
 750 c       aad(i,1)=0.3D0*4.0D0**6
 751 C "Hard" SC-p repulsion
 752         bad(i,1)=3.0D0*4.0D0**6
 753 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 754         bad(i,2)=bad(i,1)
 755 c       aad(i,1)=0.0D0
 756 c       aad(i,2)=0.0D0
 757 c       bad(i,1)=1228.8D0
 758 c       bad(i,2)=1228.8D0
 759       enddo
 760 #else
 761 C
 762 C 8/9/01 Read the SC-p interaction constants from file
 763 C
 764       do i=1,ntyp
 765         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 766       enddo
 767       do i=1,ntyp
 768         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 769         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 770         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 771         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 772       enddo
 773
 774       if (lprint) then
 775         write (iout,*) "Parameters of SC-p interactions:"
 776         do i=1,20
 777           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 778      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 779         enddo
 780       endif
 781 #endif
 782 C
 783 C Define the constants of the disulfide bridge
 784 C
 785       ebr=-5.50D0
 786 c
 787 c Old arbitrary potential - commented out.
 788 c
 789 c      dbr= 4.20D0
 790 c      fbr= 3.30D0
 791 c
 792 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 793 c energy surface of diethyl disulfide.
 794 c A. Liwo and U. Kozlowska, 11/24/03
 795 c
 796       D0CM = 3.78d0
 797       AKCM = 15.1d0
 798       AKTH = 11.0d0
 799       AKCT = 12.0d0
 800       V1SS =-1.08d0
 801       V2SS = 7.61d0
 802       V3SS = 13.7d0
 803
 804       write (iout,'(/a)') "Disulfide bridge parameters:"
 805       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 806       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 807       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 808       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 809      & ' v3ss:',v3ss
 810       return
 811       end