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