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

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