source/cluster/wham/src-old/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 #ifdef AIX
 250       call flush_(iout)
 251 #else
 252       call flush(iout)
 253 #endif
 254       endif
 255 #endif
 256
 257 #ifdef CRYST_SC
 258 C
 259 C Read the parameters of the probability distribution/energy expression
 260 C of the side chains.
 261 C
 262       do i=1,ntyp
 263         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 264         if (i.eq.10) then
 265           dsc_inv(i)=0.0D0
 266         else
 267           dsc_inv(i)=1.0D0/dsc(i)
 268         endif
 269         if (i.ne.10) then
 270         do j=1,nlob(i)
 271           do k=1,3
 272             do l=1,3
 273               blower(l,k,j)=0.0D0
 274             enddo
 275           enddo
 276         enddo
 277         bsc(1,i)=0.0D0
 278         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 279         do j=2,nlob(i)
 280           read (irotam,*) bsc(j,i)
 281           read (irotam,*) (censc(k,j,i),k=1,3),
 282      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 283         enddo
 284         do j=1,nlob(i)
 285           do k=1,3
 286             do l=1,k
 287               akl=0.0D0
 288               do m=1,3
 289                 akl=akl+blower(k,m,j)*blower(l,m,j)
 290               enddo
 291               gaussc(k,l,j,i)=akl
 292               gaussc(l,k,j,i)=akl
 293             enddo
 294           enddo
 295         enddo
 296         endif
 297       enddo
 298       close (irotam)
 299       if (lprint) then
 300         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 301         do i=1,ntyp
 302           nlobi=nlob(i)
 303           if (nlobi.gt.0) then
 304           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 305      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 306 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 307 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 308 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 309 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 310            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 311      &                             'log h',(bsc(j,i),j=1,nlobi)
 312            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 313      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 314 c          write (iout,'(a)')
 315 c         do j=1,nlobi
 316 c           ind=0
 317 c           do k=1,3
 318 c             do l=1,k
 319 c              ind=ind+1
 320 c              blower(k,l,j)=gaussc(ind,j,i)
 321 c             enddo
 322 c           enddo
 323 c         enddo
 324           do k=1,3
 325             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 326      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 327           enddo
 328           endif
 329         enddo
 330       endif
 331 #else
 332 C
 333 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 334 C added by Urszula Kozlowska 07/11/2007
 335 C
 336       do i=1,ntyp
 337         read (irotam,*)
 338        if (i.eq.10) then
 339          read (irotam,*)
 340        else
 341          do j=1,65
 342            read(irotam,*) sc_parmin(j,i)
 343          enddo
 344        endif
 345       enddo
 346 #endif
 347       close(irotam)
 348 #ifdef CRYST_TOR
 349 C
 350 C Read torsional parameters in old format
 351 C
 352       read (itorp,*) ntortyp,nterm_old
 353       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 354       read (itorp,*) (itortyp(i),i=1,ntyp)
 355       do i=1,ntortyp
 356         do j=1,ntortyp
 357           read (itorp,'(a)')
 358           do k=1,nterm_old
 359             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 360           enddo
 361         enddo
 362       enddo
 363       close (itorp)
 364       if (lprint) then
 365         write (iout,'(/a/)') 'Torsional constants:'
 366         do i=1,ntortyp
 367           do j=1,ntortyp
 368             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 369             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 370           enddo
 371         enddo
 372       endif
 373 #else
 374 C
 375 C Read torsional parameters
 376 C
 377       read (itorp,*) ntortyp
 378       read (itorp,*) (itortyp(i),i=1,ntyp)
 379       write (iout,*) 'ntortyp',ntortyp
 380       do i=1,ntortyp
 381         do j=1,ntortyp
 382           read (itorp,*) nterm(i,j),nlor(i,j)
 383           v0ij=0.0d0
 384           si=-1.0d0
 385           do k=1,nterm(i,j)
 386             read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
 387             v0ij=v0ij+si*v1(k,i,j)
 388             si=-si
 389           enddo
 390           do k=1,nlor(i,j)
 391             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 392             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 393           enddo
 394           v0(i,j)=v0ij
 395         enddo
 396       enddo
 397       close (itorp)
 398       if (lprint) then
 399         write (iout,'(/a/)') 'Torsional constants:'
 400         do i=1,ntortyp
 401           do j=1,ntortyp
 402             write (iout,*) 'ityp',i,' jtyp',j
 403             write (iout,*) 'Fourier constants'
 404             do k=1,nterm(i,j)
 405               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 406             enddo
 407             write (iout,*) 'Lorenz constants'
 408             do k=1,nlor(i,j)
 409               write (iout,'(3(1pe15.5))')
 410      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 411             enddo
 412           enddo
 413         enddo
 414       endif
 415 C
 416 C 6/23/01 Read parameters for double torsionals
 417 C
 418       do i=1,ntortyp
 419         do j=1,ntortyp
 420           do k=1,ntortyp
 421             read (itordp,'(3a1)') t1,t2,t3
 422             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 423      &        .or. t3.ne.onelett(k)) then
 424               write (iout,*) "Error in double torsional parameter file",
 425      &         i,j,k,t1,t2,t3
 426                stop "Error in double torsional parameter file"
 427             endif
 428             read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
 429             read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
 430             read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
 431             read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
 432             read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
 433             read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
 434      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
 435           enddo
 436         enddo
 437       enddo
 438       if (lprint) then
 439       write (iout,*)
 440       write (iout,*) 'Constants for double torsionals'
 441       do i=1,ntortyp
 442         do j=1,ntortyp
 443           do k=1,ntortyp
 444             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 445      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 446             write (iout,*)
 447             write (iout,*) 'Single angles:'
 448             do l=1,ntermd_1(i,j,k)
 449               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 450      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 451      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 452             enddo
 453             write (iout,*)
 454             write (iout,*) 'Pairs of angles:'
 455             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 456             do l=1,ntermd_2(i,j,k)
 457               write (iout,'(i5,20f10.5)')
 458      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 459             enddo
 460             write (iout,*)
 461             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 462             do l=1,ntermd_2(i,j,k)
 463               write (iout,'(i5,20f10.5)')
 464      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 465             enddo
 466             write (iout,*)
 467           enddo
 468         enddo
 469       enddo
 470       endif
 471 #endif
 472 C Read of Side-chain backbone correlation parameters
 473 C Modified 11 May 2012 by Adasko
 474 CCC
 475 C
 476       read (isccor,*) nsccortyp
 477       read (isccor,*) (isccortyp(i),i=1,ntyp)
 478 c      write (iout,*) 'ntortyp',ntortyp
 479       maxinter=3
 480 cc maxinter is maximum interaction sites
 481       do l=1,maxinter
 482       do i=1,nsccortyp
 483         do j=1,nsccortyp
 484           read (isccor,*)
 485      &nterm_sccor(i,j),nlor_sccor(i,j)
 486           v0ijsccor=0.0d0
 487           si=-1.0d0
 488
 489           do k=1,nterm_sccor(i,j)
 490             read (isccor,*) kk,v1sccor(k,l,i,j)
 491      &    ,v2sccor(k,l,i,j)
 492             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 493             si=-si
 494           enddo
 495           do k=1,nlor_sccor(i,j)
 496             read (isccor,*) kk,vlor1sccor(k,i,j),
 497      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 498             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 499      &(1+vlor3sccor(k,i,j)**2)
 500           enddo
 501           v0sccor(i,j)=v0ijsccor
 502         enddo
 503       enddo
 504       enddo
 505       close (isccor)
 506
 507       if (lprint) then
 508         write (iout,'(/a/)') 'Torsional constants:'
 509         do l=1,maxinter
 510         do i=1,nsccortyp
 511           do j=1,nsccortyp
 512             write (iout,*) 'ityp',i,' jtyp',j
 513             write (iout,*) 'Fourier constants'
 514             do k=1,nterm_sccor(i,j)
 515               write (iout,'(2(1pe15.5))')
 516      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 517             enddo
 518             write (iout,*) 'Lorenz constants'
 519             do k=1,nlor_sccor(i,j)
 520               write (iout,'(3(1pe15.5))')
 521      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 522             enddo
 523           enddo
 524         enddo
 525         enddo
 526       endif
 527
 528 C
 529 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 530 C         interaction energy of the Gly, Ala, and Pro prototypes.
 531 C
 532       read (ifourier,*) nloctyp
 533       do i=1,nloctyp
 534         read (ifourier,*)
 535         read (ifourier,*) (b(ii,i),ii=1,13)
 536         if (lprint) then
 537         write (iout,*) 'Type',i
 538         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 539         endif
 540         B1(1,i)  = b(3,i)
 541         B1(2,i)  = b(5,i)
 542         B1tilde(1,i) = b(3,i)
 543         B1tilde(2,i) =-b(5,i)
 544         B2(1,i)  = b(2,i)
 545         B2(2,i)  = b(4,i)
 546         CC(1,1,i)= b(7,i)
 547         CC(2,2,i)=-b(7,i)
 548         CC(2,1,i)= b(9,i)
 549         CC(1,2,i)= b(9,i)
 550         Ctilde(1,1,i)=b(7,i)
 551         Ctilde(1,2,i)=b(9,i)
 552         Ctilde(2,1,i)=-b(9,i)
 553         Ctilde(2,2,i)=b(7,i)
 554         DD(1,1,i)= b(6,i)
 555         DD(2,2,i)=-b(6,i)
 556         DD(2,1,i)= b(8,i)
 557         DD(1,2,i)= b(8,i)
 558         Dtilde(1,1,i)=b(6,i)
 559         Dtilde(1,2,i)=b(8,i)
 560         Dtilde(2,1,i)=-b(8,i)
 561         Dtilde(2,2,i)=b(6,i)
 562         EE(1,1,i)= b(10,i)+b(11,i)
 563         EE(2,2,i)=-b(10,i)+b(11,i)
 564         EE(2,1,i)= b(12,i)-b(13,i)
 565         EE(1,2,i)= b(12,i)+b(13,i)
 566       enddo
 567       if (lprint) then
 568       do i=1,nloctyp
 569         write (iout,*) 'Type',i
 570         write (iout,*) 'B1'
 571 c        write (iout,'(f10.5)') B1(:,i)
 572         write(iout,*) B1(1,i),B1(2,i)
 573         write (iout,*) 'B2'
 574 c        write (iout,'(f10.5)') B2(:,i)
 575         write(iout,*) B2(1,i),B2(2,i)
 576         write (iout,*) 'CC'
 577         do j=1,2
 578           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 579         enddo
 580         write(iout,*) 'DD'
 581         do j=1,2
 582           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 583         enddo
 584         write(iout,*) 'EE'
 585         do j=1,2
 586           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 587         enddo
 588       enddo
 589       endif
 590 C
 591 C Read electrostatic-interaction parameters
 592 C
 593       if (lprint) then
 594         write (iout,'(/a)') 'Electrostatic interaction constants:'
 595         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 596      &            'IT','JT','APP','BPP','AEL6','AEL3'
 597       endif
 598       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 599       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 600       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 601       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 602       close (ielep)
 603       do i=1,2
 604         do j=1,2
 605         rri=rpp(i,j)**6
 606         app (i,j)=epp(i,j)*rri*rri
 607         bpp (i,j)=-2.0D0*epp(i,j)*rri
 608         ael6(i,j)=elpp6(i,j)*4.2D0**6
 609         ael3(i,j)=elpp3(i,j)*4.2D0**3
 610         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 611      &                    ael6(i,j),ael3(i,j)
 612         enddo
 613       enddo
 614 C
 615 C Read side-chain interaction parameters.
 616 C
 617       read (isidep,*) ipot,expon
 618       if (ipot.lt.1 .or. ipot.gt.5) then
 619         write (iout,'(2a)') 'Error while reading SC interaction',
 620      &               'potential file - unknown potential type.'
 621         stop
 622       endif
 623       expon2=expon/2
 624       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 625      & ', exponents are ',expon,2*expon
 626       goto (10,20,30,30,40) ipot
 627 C----------------------- LJ potential ---------------------------------
 628    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 629       if (lprint) then
 630         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 631         write (iout,'(a/)') 'The epsilon array:'
 632         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 633         write (iout,'(/a)') 'One-body parameters:'
 634         write (iout,'(a,4x,a)') 'residue','sigma'
 635         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 636       endif
 637       goto 50
 638 C----------------------- LJK potential --------------------------------
 639    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 640      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 641       if (lprint) then
 642         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 643         write (iout,'(a/)') 'The epsilon array:'
 644         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 645         write (iout,'(/a)') 'One-body parameters:'
 646         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 647         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 648      &        i=1,ntyp)
 649       endif
 650       goto 50
 651 C---------------------- GB or BP potential -----------------------------
 652    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 653      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 654      &  (alp(i),i=1,ntyp)
 655 C For the GB potential convert sigma'**2 into chi'
 656       if (ipot.eq.4) then
 657         do i=1,ntyp
 658           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 659         enddo
 660       endif
 661       if (lprint) then
 662         write (iout,'(/a/)') 'Parameters of the BP potential:'
 663         write (iout,'(a/)') 'The epsilon array:'
 664         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 665         write (iout,'(/a)') 'One-body parameters:'
 666         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 667      &       '    chip  ','    alph  '
 668         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 669      &                     chip(i),alp(i),i=1,ntyp)
 670       endif
 671       goto 50
 672 C--------------------- GBV potential -----------------------------------
 673    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 674      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 675      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 676       if (lprint) then
 677         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 678         write (iout,'(a/)') 'The epsilon array:'
 679         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 680         write (iout,'(/a)') 'One-body parameters:'
 681         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 682      &      's||/s_|_^2','    chip  ','    alph  '
 683         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 684      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 685       endif
 686    50 continue
 687       close (isidep)
 688 C-----------------------------------------------------------------------
 689 C Calculate the "working" parameters of SC interactions.
 690       do i=2,ntyp
 691         do j=1,i-1
 692           eps(i,j)=eps(j,i)
 693         enddo
 694       enddo
 695       do i=1,ntyp
 696         do j=i,ntyp
 697           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 698           sigma(j,i)=sigma(i,j)
 699           rs0(i,j)=dwa16*sigma(i,j)
 700           rs0(j,i)=rs0(i,j)
 701         enddo
 702       enddo
 703       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 704      & 'Working parameters of the SC interactions:',
 705      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 706      & '  chi1   ','   chi2   '
 707       do i=1,ntyp
 708         do j=i,ntyp
 709           epsij=eps(i,j)
 710           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 711             rrij=sigma(i,j)
 712           else
 713             rrij=rr0(i)+rr0(j)
 714           endif
 715           r0(i,j)=rrij
 716           r0(j,i)=rrij
 717           rrij=rrij**expon
 718           epsij=eps(i,j)
 719           sigeps=dsign(1.0D0,epsij)
 720           epsij=dabs(epsij)
 721           aa(i,j)=epsij*rrij*rrij
 722           bb(i,j)=-sigeps*epsij*rrij
 723           aa(j,i)=aa(i,j)
 724           bb(j,i)=bb(i,j)
 725           if (ipot.gt.2) then
 726             sigt1sq=sigma0(i)**2
 727             sigt2sq=sigma0(j)**2
 728             sigii1=sigii(i)
 729             sigii2=sigii(j)
 730             ratsig1=sigt2sq/sigt1sq
 731             ratsig2=1.0D0/ratsig1
 732             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 733             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 734             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 735           else
 736             rsum_max=sigma(i,j)
 737           endif
 738 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 739             sigmaii(i,j)=rsum_max
 740             sigmaii(j,i)=rsum_max
 741 c         else
 742 c           sigmaii(i,j)=r0(i,j)
 743 c           sigmaii(j,i)=r0(i,j)
 744 c         endif
 745 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 746           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 747             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 748             augm(i,j)=epsij*r_augm**(2*expon)
 749 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 750             augm(j,i)=augm(i,j)
 751           else
 752             augm(i,j)=0.0D0
 753             augm(j,i)=0.0D0
 754           endif
 755           if (lprint) then
 756             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 757      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 758      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 759           endif
 760         enddo
 761       enddo
 762 C
 763 C Define the SC-p interaction constants
 764 C
 765       do i=1,20
 766         do j=1,2
 767           eps_scp(i,j)=-1.5d0
 768           rscp(i,j)=4.0d0
 769         enddo
 770       enddo
 771 #ifdef OLDSCP
 772       do i=1,20
 773 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 774 C helix formation)
 775 c       aad(i,1)=0.3D0*4.0D0**12
 776 C Following line for constants currently implemented
 777 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 778         aad(i,1)=1.5D0*4.0D0**12
 779 c       aad(i,1)=0.17D0*5.6D0**12
 780         aad(i,2)=aad(i,1)
 781 C "Soft" SC-p repulsion
 782         bad(i,1)=0.0D0
 783 C Following line for constants currently implemented
 784 c       aad(i,1)=0.3D0*4.0D0**6
 785 C "Hard" SC-p repulsion
 786         bad(i,1)=3.0D0*4.0D0**6
 787 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 788         bad(i,2)=bad(i,1)
 789 c       aad(i,1)=0.0D0
 790 c       aad(i,2)=0.0D0
 791 c       bad(i,1)=1228.8D0
 792 c       bad(i,2)=1228.8D0
 793       enddo
 794 #else
 795 C
 796 C 8/9/01 Read the SC-p interaction constants from file
 797 C
 798       do i=1,ntyp
 799         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 800       enddo
 801       do i=1,ntyp
 802         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 803         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 804         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 805         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 806       enddo
 807
 808       if (lprint) then
 809         write (iout,*) "Parameters of SC-p interactions:"
 810         do i=1,20
 811           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 812      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 813         enddo
 814       endif
 815 #endif
 816 C
 817 C Define the constants of the disulfide bridge
 818 C
 819 C      ebr=-5.50D0
 820 c
 821 c Old arbitrary potential - commented out.
 822 c
 823 c      dbr= 4.20D0
 824 c      fbr= 3.30D0
 825 c
 826 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 827 c energy surface of diethyl disulfide.
 828 c A. Liwo and U. Kozlowska, 11/24/03
 829 c
 830 C      D0CM = 3.78d0
 831 C      AKCM = 15.1d0
 832 C      AKTH = 11.0d0
 833 C      AKCT = 12.0d0
 834 C      V1SS =-1.08d0
 835 C      V2SS = 7.61d0
 836 C      V3SS = 13.7d0
 837       return
 838       end