source/cluster/wham/src/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 Read of Side-chain backbone correlation parameters
 469 C Modified 11 May 2012 by Adasko
 470 CCC
 471 C
 472       read (isccor,*) nsccortyp
 473       read (isccor,*) (isccortyp(i),i=1,ntyp)
 474 c      write (iout,*) 'ntortyp',ntortyp
 475       maxinter=3
 476 cc maxinter is maximum interaction sites
 477       do l=1,maxinter
 478       do i=1,nsccortyp
 479         do j=1,nsccortyp
 480           read (isccor,*)
 481      &nterm_sccor(i,j),nlor_sccor(i,j)
 482           v0ijsccor=0.0d0
 483           si=-1.0d0
 484
 485           do k=1,nterm_sccor(i,j)
 486             read (isccor,*) kk,v1sccor(k,l,i,j)
 487      &    ,v2sccor(k,l,i,j)
 488             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 489             si=-si
 490           enddo
 491           do k=1,nlor_sccor(i,j)
 492             read (isccor,*) kk,vlor1sccor(k,i,j),
 493      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 494             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 495      &(1+vlor3sccor(k,i,j)**2)
 496           enddo
 497           v0sccor(i,j)=v0ijsccor
 498         enddo
 499       enddo
 500       enddo
 501       close (isccor)
 502
 503       if (lprint) then
 504         write (iout,'(/a/)') 'Torsional constants:'
 505         do l=1,maxinter
 506         do i=1,nsccortyp
 507           do j=1,nsccortyp
 508             write (iout,*) 'ityp',i,' jtyp',j
 509             write (iout,*) 'Fourier constants'
 510             do k=1,nterm_sccor(i,j)
 511               write (iout,'(2(1pe15.5))')
 512      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 513             enddo
 514             write (iout,*) 'Lorenz constants'
 515             do k=1,nlor_sccor(i,j)
 516               write (iout,'(3(1pe15.5))')
 517      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 518             enddo
 519           enddo
 520         enddo
 521         enddo
 522       endif
 523
 524 C
 525 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 526 C         interaction energy of the Gly, Ala, and Pro prototypes.
 527 C
 528       read (ifourier,*) nloctyp
 529       do i=1,nloctyp
 530         read (ifourier,*)
 531         read (ifourier,*) (b(ii,i),ii=1,13)
 532         if (lprint) then
 533         write (iout,*) 'Type',i
 534         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 535         endif
 536         B1(1,i)  = b(3,i)
 537         B1(2,i)  = b(5,i)
 538         B1tilde(1,i) = b(3,i)
 539         B1tilde(2,i) =-b(5,i)
 540         B2(1,i)  = b(2,i)
 541         B2(2,i)  = b(4,i)
 542         CC(1,1,i)= b(7,i)
 543         CC(2,2,i)=-b(7,i)
 544         CC(2,1,i)= b(9,i)
 545         CC(1,2,i)= b(9,i)
 546         Ctilde(1,1,i)=b(7,i)
 547         Ctilde(1,2,i)=b(9,i)
 548         Ctilde(2,1,i)=-b(9,i)
 549         Ctilde(2,2,i)=b(7,i)
 550         DD(1,1,i)= b(6,i)
 551         DD(2,2,i)=-b(6,i)
 552         DD(2,1,i)= b(8,i)
 553         DD(1,2,i)= b(8,i)
 554         Dtilde(1,1,i)=b(6,i)
 555         Dtilde(1,2,i)=b(8,i)
 556         Dtilde(2,1,i)=-b(8,i)
 557         Dtilde(2,2,i)=b(6,i)
 558         EE(1,1,i)= b(10,i)+b(11,i)
 559         EE(2,2,i)=-b(10,i)+b(11,i)
 560         EE(2,1,i)= b(12,i)-b(13,i)
 561         EE(1,2,i)= b(12,i)+b(13,i)
 562       enddo
 563       if (lprint) then
 564       do i=1,nloctyp
 565         write (iout,*) 'Type',i
 566         write (iout,*) 'B1'
 567 c        write (iout,'(f10.5)') B1(:,i)
 568         write(iout,*) B1(1,i),B1(2,i)
 569         write (iout,*) 'B2'
 570 c        write (iout,'(f10.5)') B2(:,i)
 571         write(iout,*) B2(1,i),B2(2,i)
 572         write (iout,*) 'CC'
 573         do j=1,2
 574           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 575         enddo
 576         write(iout,*) 'DD'
 577         do j=1,2
 578           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 579         enddo
 580         write(iout,*) 'EE'
 581         do j=1,2
 582           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 583         enddo
 584       enddo
 585       endif
 586 C
 587 C Read electrostatic-interaction parameters
 588 C
 589       if (lprint) then
 590         write (iout,'(/a)') 'Electrostatic interaction constants:'
 591         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 592      &            'IT','JT','APP','BPP','AEL6','AEL3'
 593       endif
 594       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 595       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 596       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 597       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 598       close (ielep)
 599       do i=1,2
 600         do j=1,2
 601         rri=rpp(i,j)**6
 602         app (i,j)=epp(i,j)*rri*rri
 603         bpp (i,j)=-2.0D0*epp(i,j)*rri
 604         ael6(i,j)=elpp6(i,j)*4.2D0**6
 605         ael3(i,j)=elpp3(i,j)*4.2D0**3
 606         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 607      &                    ael6(i,j),ael3(i,j)
 608         enddo
 609       enddo
 610 C
 611 C Read side-chain interaction parameters.
 612 C
 613       read (isidep,*) ipot,expon
 614       if (ipot.lt.1 .or. ipot.gt.5) then
 615         write (iout,'(2a)') 'Error while reading SC interaction',
 616      &               'potential file - unknown potential type.'
 617         stop
 618       endif
 619       expon2=expon/2
 620       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 621      & ', exponents are ',expon,2*expon
 622       goto (10,20,30,30,40) ipot
 623 C----------------------- LJ potential ---------------------------------
 624    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 625       if (lprint) then
 626         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 627         write (iout,'(a/)') 'The epsilon array:'
 628         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 629         write (iout,'(/a)') 'One-body parameters:'
 630         write (iout,'(a,4x,a)') 'residue','sigma'
 631         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 632       endif
 633       goto 50
 634 C----------------------- LJK potential --------------------------------
 635    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 636      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 637       if (lprint) then
 638         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 639         write (iout,'(a/)') 'The epsilon array:'
 640         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 641         write (iout,'(/a)') 'One-body parameters:'
 642         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 643         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 644      &        i=1,ntyp)
 645       endif
 646       goto 50
 647 C---------------------- GB or BP potential -----------------------------
 648    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 649      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 650      &  (alp(i),i=1,ntyp)
 651 C For the GB potential convert sigma'**2 into chi'
 652       if (ipot.eq.4) then
 653         do i=1,ntyp
 654           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 655         enddo
 656       endif
 657       if (lprint) then
 658         write (iout,'(/a/)') 'Parameters of the BP potential:'
 659         write (iout,'(a/)') 'The epsilon array:'
 660         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 661         write (iout,'(/a)') 'One-body parameters:'
 662         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 663      &       '    chip  ','    alph  '
 664         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 665      &                     chip(i),alp(i),i=1,ntyp)
 666       endif
 667       goto 50
 668 C--------------------- GBV potential -----------------------------------
 669    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 670      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 671      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 672       if (lprint) then
 673         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 674         write (iout,'(a/)') 'The epsilon array:'
 675         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 676         write (iout,'(/a)') 'One-body parameters:'
 677         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 678      &      's||/s_|_^2','    chip  ','    alph  '
 679         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 680      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 681       endif
 682    50 continue
 683       close (isidep)
 684 C-----------------------------------------------------------------------
 685 C Calculate the "working" parameters of SC interactions.
 686       do i=2,ntyp
 687         do j=1,i-1
 688           eps(i,j)=eps(j,i)
 689         enddo
 690       enddo
 691       do i=1,ntyp
 692         do j=i,ntyp
 693           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 694           sigma(j,i)=sigma(i,j)
 695           rs0(i,j)=dwa16*sigma(i,j)
 696           rs0(j,i)=rs0(i,j)
 697         enddo
 698       enddo
 699       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 700      & 'Working parameters of the SC interactions:',
 701      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 702      & '  chi1   ','   chi2   '
 703       do i=1,ntyp
 704         do j=i,ntyp
 705           epsij=eps(i,j)
 706           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 707             rrij=sigma(i,j)
 708           else
 709             rrij=rr0(i)+rr0(j)
 710           endif
 711           r0(i,j)=rrij
 712           r0(j,i)=rrij
 713           rrij=rrij**expon
 714           epsij=eps(i,j)
 715           sigeps=dsign(1.0D0,epsij)
 716           epsij=dabs(epsij)
 717           aa(i,j)=epsij*rrij*rrij
 718           bb(i,j)=-sigeps*epsij*rrij
 719           aa(j,i)=aa(i,j)
 720           bb(j,i)=bb(i,j)
 721           if (ipot.gt.2) then
 722             sigt1sq=sigma0(i)**2
 723             sigt2sq=sigma0(j)**2
 724             sigii1=sigii(i)
 725             sigii2=sigii(j)
 726             ratsig1=sigt2sq/sigt1sq
 727             ratsig2=1.0D0/ratsig1
 728             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 729             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 730             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 731           else
 732             rsum_max=sigma(i,j)
 733           endif
 734 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 735             sigmaii(i,j)=rsum_max
 736             sigmaii(j,i)=rsum_max
 737 c         else
 738 c           sigmaii(i,j)=r0(i,j)
 739 c           sigmaii(j,i)=r0(i,j)
 740 c         endif
 741 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 742           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 743             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 744             augm(i,j)=epsij*r_augm**(2*expon)
 745 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 746             augm(j,i)=augm(i,j)
 747           else
 748             augm(i,j)=0.0D0
 749             augm(j,i)=0.0D0
 750           endif
 751           if (lprint) then
 752             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 753      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 754      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 755           endif
 756         enddo
 757       enddo
 758 C
 759 C Define the SC-p interaction constants
 760 C
 761       do i=1,20
 762         do j=1,2
 763           eps_scp(i,j)=-1.5d0
 764           rscp(i,j)=4.0d0
 765         enddo
 766       enddo
 767 #ifdef OLDSCP
 768       do i=1,20
 769 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 770 C helix formation)
 771 c       aad(i,1)=0.3D0*4.0D0**12
 772 C Following line for constants currently implemented
 773 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 774         aad(i,1)=1.5D0*4.0D0**12
 775 c       aad(i,1)=0.17D0*5.6D0**12
 776         aad(i,2)=aad(i,1)
 777 C "Soft" SC-p repulsion
 778         bad(i,1)=0.0D0
 779 C Following line for constants currently implemented
 780 c       aad(i,1)=0.3D0*4.0D0**6
 781 C "Hard" SC-p repulsion
 782         bad(i,1)=3.0D0*4.0D0**6
 783 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 784         bad(i,2)=bad(i,1)
 785 c       aad(i,1)=0.0D0
 786 c       aad(i,2)=0.0D0
 787 c       bad(i,1)=1228.8D0
 788 c       bad(i,2)=1228.8D0
 789       enddo
 790 #else
 791 C
 792 C 8/9/01 Read the SC-p interaction constants from file
 793 C
 794       do i=1,ntyp
 795         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 796       enddo
 797       do i=1,ntyp
 798         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 799         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 800         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 801         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 802       enddo
 803
 804       if (lprint) then
 805         write (iout,*) "Parameters of SC-p interactions:"
 806         do i=1,20
 807           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 808      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 809         enddo
 810       endif
 811 #endif
 812 C
 813 C Define the constants of the disulfide bridge
 814 C
 815 C      ebr=-5.50D0
 816 c
 817 c Old arbitrary potential - commented out.
 818 c
 819 c      dbr= 4.20D0
 820 c      fbr= 3.30D0
 821 c
 822 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 823 c energy surface of diethyl disulfide.
 824 c A. Liwo and U. Kozlowska, 11/24/03
 825 c
 826 C      D0CM = 3.78d0
 827 C      AKCM = 15.1d0
 828 C      AKTH = 11.0d0
 829 C      AKCT = 12.0d0
 830 C      V1SS =-1.08d0
 831 C      V2SS = 7.61d0
 832 C      V3SS = 13.7d0
 833       return
 834       end