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