source/cluster/wham/src-M-NEWCORR/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       character*1 toronelet(-2:2) /"p","a","G","A","P"/
  23       logical lprint
  24       dimension blower(3,3,maxlob)
  25       double precision ip,mp
  26 C
  27 C Body
  28 C
  29 C Set LPRINT=.TRUE. for debugging
  30       dwa16=2.0d0**(1.0d0/6.0d0)
  31       lprint=.false.
  32       itypro=20
  33 C Assign virtual-bond length
  34       vbl=3.8D0
  35       vblinv=1.0D0/vbl
  36       vblinv2=vblinv*vblinv
  37 c
  38 c Read the virtual-bond parameters, masses, and moments of inertia
  39 c and Stokes' radii of the peptide group and side chains
  40 c
  41 #ifdef CRYST_BOND
  42       read (ibond,*) vbldp0,akp
  43       do i=1,ntyp
  44         nbondterm(i)=1
  45         read (ibond,*) vbldsc0(1,i),aksc(1,i)
  46         dsc(i) = vbldsc0(1,i)
  47         if (i.eq.10) then
  48           dsc_inv(i)=0.0D0
  49         else
  50           dsc_inv(i)=1.0D0/dsc(i)
  51         endif
  52       enddo
  53 #else
  54       read (ibond,*) ijunk,vbldp0,akp,rjunk
  55       do i=1,ntyp
  56         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
  57      &   j=1,nbondterm(i))
  58         dsc(i) = vbldsc0(1,i)
  59         if (i.eq.10) then
  60           dsc_inv(i)=0.0D0
  61         else
  62           dsc_inv(i)=1.0D0/dsc(i)
  63         endif
  64       enddo
  65 #endif
  66       if (lprint) then
  67         write(iout,'(/a/)')"Force constants virtual bonds:"
  68         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
  69      &   'inertia','Pstok'
  70         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
  71         do i=1,ntyp
  72           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
  73      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
  74           do j=2,nbondterm(i)
  75             write (iout,'(13x,3f10.5)')
  76      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
  77           enddo
  78         enddo
  79       endif
  80 #ifdef CRYST_THETA
  81 C
  82 C Read the parameters of the probability distribution/energy expression
  83 C of the virtual-bond valence angles theta
  84 C
  85       do i=1,ntyp
  86         read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
  87      &    (bthet(j,i,1,1),j=1,2)
  88         read (ithep,*) (polthet(j,i),j=0,3)
  89         read (ithep,*) (gthet(j,i),j=1,3)
  90         read (ithep,*) theta0(i),sig0(i),sigc0(i)
  91         sigc0(i)=sigc0(i)**2
  92       enddo
  93       do i=1,ntyp
  94       athet(1,i,1,-1)=athet(1,i,1,1)
  95       athet(2,i,1,-1)=athet(2,i,1,1)
  96       bthet(1,i,1,-1)=-bthet(1,i,1,1)
  97       bthet(2,i,1,-1)=-bthet(2,i,1,1)
  98       athet(1,i,-1,1)=-athet(1,i,1,1)
  99       athet(2,i,-1,1)=-athet(2,i,1,1)
 100       bthet(1,i,-1,1)=bthet(1,i,1,1)
 101       bthet(2,i,-1,1)=bthet(2,i,1,1)
 102       enddo
 103       do i=-ntyp,-1
 104       a0thet(i)=a0thet(-i)
 105       athet(1,i,-1,-1)=athet(1,-i,1,1)
 106       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 107       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 108       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 109       athet(1,i,-1,1)=athet(1,-i,1,1)
 110       athet(2,i,-1,1)=-athet(2,-i,1,1)
 111       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 112       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 113       athet(1,i,1,-1)=-athet(1,-i,1,1)
 114       athet(2,i,1,-1)=athet(2,-i,1,1)
 115       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 116       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 117       theta0(i)=theta0(-i)
 118       sig0(i)=sig0(-i)
 119       sigc0(i)=sigc0(-i)
 120        do j=0,3
 121         polthet(j,i)=polthet(j,-i)
 122        enddo
 123        do j=1,3
 124          gthet(j,i)=gthet(j,-i)
 125        enddo
 126       enddo
 127       close (ithep)
 128       if (lprint) then
 129 c       write (iout,'(a)')
 130 c    &    'Parameters of the virtual-bond valence angles:'
 131 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 132 c    & '    ATHETA0   ','         A1   ','        A2    ',
 133 c    & '        B1    ','         B2   '
 134 c       do i=1,ntyp
 135 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 136 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 137 c       enddo
 138 c       write (iout,'(/a/9x,5a/79(1h-))')
 139 c    & 'Parameters of the expression for sigma(theta_c):',
 140 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 141 c    & '     ALPH3    ','    SIGMA0C   '
 142 c       do i=1,ntyp
 143 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 144 c    &      (polthet(j,i),j=0,3),sigc0(i)
 145 c       enddo
 146 c       write (iout,'(/a/9x,5a/79(1h-))')
 147 c    & 'Parameters of the second gaussian:',
 148 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 149 c    & '        G2    ','         G3   '
 150 c       do i=1,ntyp
 151 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 152 c    &       sig0(i),(gthet(j,i),j=1,3)
 153 c       enddo
 154         write (iout,'(a)')
 155      &    'Parameters of the virtual-bond valence angles:'
 156         write (iout,'(/a/9x,5a/79(1h-))')
 157      & 'Coefficients of expansion',
 158      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 159      & '   b1*10^1    ','    b2*10^1   '
 160         do i=1,ntyp
 161           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 162      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
 163      &        (10*bthet(j,i,1,1),j=1,2)
 164         enddo
 165         write (iout,'(/a/9x,5a/79(1h-))')
 166      & 'Parameters of the expression for sigma(theta_c):',
 167      & ' alpha0       ','  alph1       ',' alph2        ',
 168      & ' alhp3        ','   sigma0c    '
 169         do i=1,ntyp
 170           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 171      &      (polthet(j,i),j=0,3),sigc0(i)
 172         enddo
 173         write (iout,'(/a/9x,5a/79(1h-))')
 174      & 'Parameters of the second gaussian:',
 175      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 176      & '        G2    ','   G3*10^1    '
 177         do i=1,ntyp
 178           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 179      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 180         enddo
 181       endif
 182 #else
 183 C
 184 C Read the parameters of Utheta determined from ab initio surfaces
 185 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 186 C
 187 c      write (iout,*) "tu dochodze"
 188       read (ithep,*) nthetyp,ntheterm,ntheterm2,
 189      &  ntheterm3,nsingle,ndouble
 190       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 191       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 192       do i=-ntyp1,-1
 193         ithetyp(i)=-ithetyp(-i)
 194       enddo
 195 c      write (iout,*) "tu dochodze"
 196       do iblock=1,2
 197       do i=-maxthetyp,maxthetyp
 198         do j=-maxthetyp,maxthetyp
 199           do k=-maxthetyp,maxthetyp
 200             aa0thet(i,j,k,iblock)=0.0d0
 201             do l=1,ntheterm
 202               aathet(l,i,j,k,iblock)=0.0d0
 203             enddo
 204             do l=1,ntheterm2
 205               do m=1,nsingle
 206                 bbthet(m,l,i,j,k,iblock)=0.0d0
 207                 ccthet(m,l,i,j,k,iblock)=0.0d0
 208                 ddthet(m,l,i,j,k,iblock)=0.0d0
 209                 eethet(m,l,i,j,k,iblock)=0.0d0
 210               enddo
 211             enddo
 212             do l=1,ntheterm3
 213               do m=1,ndouble
 214                 do mm=1,ndouble
 215                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 216                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 217                 enddo
 218               enddo
 219             enddo
 220           enddo
 221         enddo
 222       enddo
 223       enddo
 224       do iblock=1,2
 225       do i=0,nthetyp
 226         do j=-nthetyp,nthetyp
 227           do k=-nthetyp,nthetyp
 228             read (ithep,'(6a)') res1
 229             read (ithep,*) aa0thet(i,j,k,iblock)
 230             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 231             read (ithep,*)
 232      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 233      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 234      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 235      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 236      &        ,ll=1,ntheterm2)
 237             read (ithep,*)
 238      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 239      &      ffthet(lll,llll,ll,i,j,k,iblock),
 240      &         ggthet(llll,lll,ll,i,j,k,iblock)
 241      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 242      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 243           enddo
 244         enddo
 245       enddo
 246 C
 247 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 248 C coefficients of theta-and-gamma-dependent terms are zero.
 249 C
 250       do i=1,nthetyp
 251         do j=1,nthetyp
 252           do l=1,ntheterm
 253             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 254             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 255           enddo
 256           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 257           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 258         enddo
 259         do l=1,ntheterm
 260           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 261         enddo
 262         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 263       enddo
 264       enddo
 265 C Substitution for D aminoacids from symmetry.
 266       do iblock=1,2
 267       do i=-nthetyp,0
 268         do j=-nthetyp,nthetyp
 269           do k=-nthetyp,nthetyp
 270            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 271            do l=1,ntheterm
 272            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 273            enddo
 274            do ll=1,ntheterm2
 275             do lll=1,nsingle
 276             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 277             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 278             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 279             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 280             enddo
 281           enddo
 282           do ll=1,ntheterm3
 283            do lll=2,ndouble
 284             do llll=1,lll-1
 285             ffthet(llll,lll,ll,i,j,k,iblock)=
 286      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 287             ffthet(lll,llll,ll,i,j,k,iblock)=
 288      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 289             ggthet(llll,lll,ll,i,j,k,iblock)=
 290      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 291             ggthet(lll,llll,ll,i,j,k,iblock)=
 292      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 293             enddo !ll
 294            enddo  !lll
 295           enddo   !llll
 296          enddo    !k
 297         enddo     !j
 298        enddo      !i
 299       enddo       !iblock
 300
 301 C
 302 C Control printout of the coefficients of virtual-bond-angle potentials
 303 C
 304       if (lprint) then
 305         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 306         do iblock=1,2
 307         do i=1,nthetyp+1
 308           do j=1,nthetyp+1
 309             do k=1,nthetyp+1
 310               write (iout,'(//4a)')
 311      &         'Type ',onelett(i),onelett(j),onelett(k)
 312               write (iout,'(//a,10x,a)') " l","a[l]"
 313               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 314               write (iout,'(i2,1pe15.5)')
 315      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 316             do l=1,ntheterm2
 317               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 318      &          "b",l,"c",l,"d",l,"e",l
 319               do m=1,nsingle
 320                 write (iout,'(i2,4(1pe15.5))') m,
 321      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 322      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 323               enddo
 324             enddo
 325             do l=1,ntheterm3
 326               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 327      &          "f+",l,"f-",l,"g+",l,"g-",l
 328               do m=2,ndouble
 329                 do n=1,m-1
 330                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 331      &              ffthet(n,m,l,i,j,k,iblock),
 332      &              ffthet(m,n,l,i,j,k,iblock),
 333      &              ggthet(n,m,l,i,j,k,iblock),
 334      &              ggthet(m,n,l,i,j,k,iblock)
 335                 enddo
 336               enddo
 337             enddo
 338           enddo
 339         enddo
 340         enddo
 341       enddo
 342       call flush(iout)
 343       endif
 344 #endif
 345
 346 #ifdef CRYST_SC
 347 C
 348 C Read the parameters of the probability distribution/energy expression
 349 C of the side chains.
 350 C
 351       do i=1,ntyp
 352 cc      write (iout,*) "tu dochodze",i
 353         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 354         if (i.eq.10) then
 355           dsc_inv(i)=0.0D0
 356         else
 357           dsc_inv(i)=1.0D0/dsc(i)
 358         endif
 359         if (i.ne.10) then
 360         do j=1,nlob(i)
 361           do k=1,3
 362             do l=1,3
 363               blower(l,k,j)=0.0D0
 364             enddo
 365           enddo
 366         enddo
 367         bsc(1,i)=0.0D0
 368         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 369         censc(1,1,-i)=censc(1,1,i)
 370         censc(2,1,-i)=censc(2,1,i)
 371         censc(3,1,-i)=-censc(3,1,i)
 372         do j=2,nlob(i)
 373           read (irotam,*) bsc(j,i)
 374           read (irotam,*) (censc(k,j,i),k=1,3),
 375      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 376         censc(1,j,-i)=censc(1,j,i)
 377         censc(2,j,-i)=censc(2,j,i)
 378         censc(3,j,-i)=-censc(3,j,i)
 379 C BSC is amplitude of Gaussian
 380         enddo
 381         do j=1,nlob(i)
 382           do k=1,3
 383             do l=1,k
 384               akl=0.0D0
 385               do m=1,3
 386                 akl=akl+blower(k,m,j)*blower(l,m,j)
 387               enddo
 388               gaussc(k,l,j,i)=akl
 389               gaussc(l,k,j,i)=akl
 390              if (((k.eq.3).and.(l.ne.3))
 391      &        .or.((l.eq.3).and.(k.ne.3))) then
 392                 gaussc(k,l,j,-i)=-akl
 393                 gaussc(l,k,j,-i)=-akl
 394               else
 395                 gaussc(k,l,j,-i)=akl
 396                 gaussc(l,k,j,-i)=akl
 397               endif
 398             enddo
 399           enddo
 400         enddo
 401         endif
 402       enddo
 403       close (irotam)
 404       if (lprint) then
 405         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 406         do i=1,ntyp
 407           nlobi=nlob(i)
 408           if (nlobi.gt.0) then
 409           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 410      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 411 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 412 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 413 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 414 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 415            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 416      &                             'log h',(bsc(j,i),j=1,nlobi)
 417            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 418      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 419 c          write (iout,'(a)')
 420 c         do j=1,nlobi
 421 c           ind=0
 422 c           do k=1,3
 423 c             do l=1,k
 424 c              ind=ind+1
 425 c              blower(k,l,j)=gaussc(ind,j,i)
 426 c             enddo
 427 c           enddo
 428 c         enddo
 429           do k=1,3
 430             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 431      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 432           enddo
 433           endif
 434         enddo
 435       endif
 436 #else
 437 C
 438 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 439 C added by Urszula Kozlowska 07/11/2007
 440 C
 441       do i=1,ntyp
 442         read (irotam,*)
 443        if (i.eq.10) then
 444          read (irotam,*)
 445        else
 446          do j=1,65
 447            read(irotam,*) sc_parmin(j,i)
 448          enddo
 449        endif
 450       enddo
 451 #endif
 452       close(irotam)
 453 #ifdef CRYST_TOR
 454 C
 455 C Read torsional parameters in old format
 456 C
 457       read (itorp,*) ntortyp,nterm_old
 458       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 459       read (itorp,*) (itortyp(i),i=1,ntyp)
 460       do i=1,ntortyp
 461         do j=1,ntortyp
 462           read (itorp,'(a)')
 463           do k=1,nterm_old
 464             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 465           enddo
 466         enddo
 467       enddo
 468       close (itorp)
 469       if (lprint) then
 470         write (iout,'(/a/)') 'Torsional constants:'
 471         do i=1,ntortyp
 472           do j=1,ntortyp
 473             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 474             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 475           enddo
 476         enddo
 477       endif
 478
 479
 480 #else
 481 C
 482 C Read torsional parameters
 483 C
 484       read (itorp,*) ntortyp
 485       read (itorp,*) (itortyp(i),i=1,ntyp)
 486       write (iout,*) 'ntortyp',ntortyp
 487       do iblock=1,2
 488       do i=-ntyp,-1
 489        itortyp(i)=-itortyp(-i)
 490       enddo
 491 c      write (iout,*) 'ntortyp',ntortyp
 492       do i=0,ntortyp-1
 493         do j=-ntortyp+1,ntortyp-1
 494           read (itorp,*) nterm(i,j,iblock),
 495      &          nlor(i,j,iblock)
 496           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 497           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 498           v0ij=0.0d0
 499           si=-1.0d0
 500           do k=1,nterm(i,j,iblock)
 501             read (itorp,*) kk,v1(k,i,j,iblock),
 502      &      v2(k,i,j,iblock)
 503             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 504             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 505             v0ij=v0ij+si*v1(k,i,j,iblock)
 506             si=-si
 507          enddo
 508           do k=1,nlor(i,j,iblock)
 509             read (itorp,*) kk,vlor1(k,i,j),
 510      &        vlor2(k,i,j),vlor3(k,i,j)
 511             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 512           enddo
 513           v0(i,j,iblock)=v0ij
 514           v0(-i,-j,iblock)=v0ij
 515         enddo
 516       enddo
 517       enddo
 518       close (itorp)
 519       if (lprint) then
 520         write (iout,'(/a/)') 'Torsional constants:'
 521         do i=1,ntortyp
 522           do j=1,ntortyp
 523             write (iout,*) 'ityp',i,' jtyp',j
 524             write (iout,*) 'Fourier constants'
 525             do iblock=1,2
 526             do k=1,nterm(i,j,iblock)
 527               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 528      &        v2(k,i,j,iblock)
 529             enddo
 530             write (iout,*) 'Lorenz constants'
 531             do k=1,nlor(i,j,iblock)
 532               write (iout,'(3(1pe15.5))')
 533      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 534             enddo
 535             enddo
 536           enddo
 537         enddo
 538       endif
 539 C
 540 C 6/23/01 Read parameters for double torsionals
 541 C
 542       do iblock=1,2
 543       do i=0,ntortyp-1
 544         do j=-ntortyp+1,ntortyp-1
 545           do k=-ntortyp+1,ntortyp-1
 546             read (itordp,'(3a1)') t1,t2,t3
 547 c              write (iout,*) "OK onelett",
 548 c     &         i,j,k,t1,t2,t3
 549
 550             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 551      &        .or. t3.ne.toronelet(k)) then
 552               write (iout,*) "Error in double torsional parameter file",
 553      &         i,j,k,t1,t2,t3
 554 #ifdef MPI
 555               call MPI_Finalize(Ierror)
 556 #endif
 557                stop "Error in double torsional parameter file"
 558             endif
 559           read (itordp,*) ntermd_1(i,j,k,iblock),
 560      &         ntermd_2(i,j,k,iblock)
 561             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 562             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 563             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
 564      &         ntermd_1(i,j,k,iblock))
 565             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
 566      &         ntermd_1(i,j,k,iblock))
 567             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
 568      &         ntermd_1(i,j,k,iblock))
 569             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
 570      &         ntermd_1(i,j,k,iblock))
 571 C Martix of D parameters for one dimesional foureir series
 572             do l=1,ntermd_1(i,j,k,iblock)
 573              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 574              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 575              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 576              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 577 c            write(iout,*) "whcodze" ,
 578 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 579             enddo
 580             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
 581      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 582      &         v2s(m,l,i,j,k,iblock),
 583      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 584 C Martix of D parameters for two dimesional fourier series
 585             do l=1,ntermd_2(i,j,k,iblock)
 586              do m=1,l-1
 587              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 588              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 589              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 590              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 591              enddo!m
 592             enddo!l
 593           enddo!k
 594         enddo!j
 595       enddo!i
 596       enddo!iblock
 597       if (lprint) then
 598       write (iout,*)
 599       write (iout,*) 'Constants for double torsionals'
 600       do iblock=1,2
 601       do i=0,ntortyp-1
 602         do j=-ntortyp+1,ntortyp-1
 603           do k=-ntortyp+1,ntortyp-1
 604             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 605      &        ' nsingle',ntermd_1(i,j,k,iblock),
 606      &        ' ndouble',ntermd_2(i,j,k,iblock)
 607             write (iout,*)
 608             write (iout,*) 'Single angles:'
 609             do l=1,ntermd_1(i,j,k,iblock)
 610               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 611      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 612      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 613      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 614             enddo
 615             write (iout,*)
 616             write (iout,*) 'Pairs of angles:'
 617             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 618             do l=1,ntermd_2(i,j,k,iblock)
 619               write (iout,'(i5,20f10.5)')
 620      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 621             enddo
 622             write (iout,*)
 623            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 624             do l=1,ntermd_2(i,j,k,iblock)
 625               write (iout,'(i5,20f10.5)')
 626      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 627      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 628             enddo
 629             write (iout,*)
 630           enddo
 631         enddo
 632       enddo
 633       enddo
 634       endif
 635 #endif
 636 C Read of Side-chain backbone correlation parameters
 637 C Modified 11 May 2012 by Adasko
 638 CCC
 639 C
 640       read (isccor,*) nsccortyp
 641       read (isccor,*) (isccortyp(i),i=1,ntyp)
 642       do i=-ntyp,-1
 643         isccortyp(i)=-isccortyp(-i)
 644       enddo
 645       iscprol=isccortyp(20)
 646 c      write (iout,*) 'ntortyp',ntortyp
 647       maxinter=3
 648 cc maxinter is maximum interaction sites
 649       do l=1,maxinter
 650       do i=1,nsccortyp
 651         do j=1,nsccortyp
 652           read (isccor,*)
 653      &nterm_sccor(i,j),nlor_sccor(i,j)
 654           write (iout,*) nterm_sccor(i,j)
 655           v0ijsccor=0.0d0
 656           v0ijsccor1=0.0d0
 657           v0ijsccor2=0.0d0
 658           v0ijsccor3=0.0d0
 659           si=-1.0d0
 660           nterm_sccor(-i,j)=nterm_sccor(i,j)
 661           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 662           nterm_sccor(i,-j)=nterm_sccor(i,j)
 663           write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
 664      &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
 665           do k=1,nterm_sccor(i,j)
 666             read (isccor,*) kk,v1sccor(k,l,i,j)
 667      &    ,v2sccor(k,l,i,j)
 668             if (j.eq.iscprol) then
 669              if (i.eq.isccortyp(10)) then
 670              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 671              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 672              else
 673              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 674      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 675              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 676      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 677              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 678              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 679              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 680              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 681              endif
 682             else
 683              if (i.eq.isccortyp(10)) then
 684              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 685              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 686              else
 687                if (j.eq.isccortyp(10)) then
 688              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 689              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 690                else
 691              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 692              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 693              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 694              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 695              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 696              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 697                 endif
 698                endif
 699             endif
 700             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 701             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 702             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 703             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 704             si=-si
 705            enddo
 706           do k=1,nlor_sccor(i,j)
 707             read (isccor,*) kk,vlor1sccor(k,i,j),
 708      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 709             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 710      &(1+vlor3sccor(k,i,j)**2)
 711           enddo
 712           v0sccor(l,i,j)=v0ijsccor
 713           v0sccor(l,-i,j)=v0ijsccor1
 714           v0sccor(l,i,-j)=v0ijsccor2
 715           v0sccor(l,-i,-j)=v0ijsccor3
 716           enddo
 717         enddo
 718       enddo
 719       close (isccor)
 720       if (lprint) then
 721         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 722         do i=1,nsccortyp
 723           do j=1,nsccortyp
 724             write (iout,*) 'ityp',i,' jtyp',j
 725             write (iout,*) 'Fourier constants'
 726             do l=1,3
 727             do k=1,nterm_sccor(i,j)
 728               write (iout,'(2(1pe15.5))')
 729      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 730             enddo
 731             enddo
 732             write (iout,*) 'Lorenz constants'
 733             do k=1,nlor_sccor(i,j)
 734               write (iout,'(3(1pe15.5))')
 735      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 736             enddo
 737           enddo
 738         enddo
 739       endif
 740 C
 741 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 742 C         interaction energy of the Gly, Ala, and Pro prototypes.
 743 C
 744       write (iout,*) "Reading Fourier"
 745       call flush(iout)
 746       read (ifourier,*) nloctyp
 747       do i=0,nloctyp-1
 748         read (ifourier,*)
 749         read (ifourier,*) (b(ii,i),ii=1,13)
 750         if (lprint) then
 751         write (iout,*) 'Type',i
 752         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 753         endif
 754       write (iout,*) "Reading Fourier 1"
 755       call flush(iout)
 756 #ifdef NEWCORR
 757         read (ifourier,*) (bnew1(ii,1,i),ii=1,3)
 758         read (ifourier,*) (bnew2(ii,1,i),ii=1,3)
 759         read (ifourier,*) (bnew1(ii,2,i),ii=1,1)
 760         read (ifourier,*) (bnew2(ii,2,i),ii=1,1)
 761         read (ifourier,*) (eenew(ii,i),ii=1,1)
 762         CC(1,1,i)= b(7,i)
 763         CC(2,2,i)=-b(7,i)
 764         CC(2,1,i)= b(9,i)
 765         CC(1,2,i)= b(9,i)
 766         CC(1,1,-i)= b(7,i)
 767         CC(2,2,-i)=-b(7,i)
 768         CC(2,1,-i)=-b(9,i)
 769         CC(1,2,-i)=-b(9,i)
 770         Ctilde(1,1,i)=b(7,i)
 771         Ctilde(1,2,i)=b(9,i)
 772         Ctilde(2,1,i)=-b(9,i)
 773         Ctilde(2,2,i)=b(7,i)
 774         Ctilde(1,1,-i)=b(7,i)
 775         Ctilde(1,2,-i)=-b(9,i)
 776         Ctilde(2,1,-i)=b(9,i)
 777         Ctilde(2,2,-i)=b(7,i)
 778
 779         DD(1,1,i)= b(6,i)
 780         DD(2,2,i)=-b(6,i)
 781         DD(2,1,i)= b(8,i)
 782         DD(1,2,i)= b(8,i)
 783         DD(1,1,-i)= b(6,i)
 784         DD(2,2,-i)=-b(6,i)
 785         DD(2,1,-i)=-b(8,i)
 786         DD(1,2,-i)=-b(8,i)
 787         Dtilde(1,1,i)=b(6,i)
 788         Dtilde(1,2,i)=b(8,i)
 789         Dtilde(2,1,i)=-b(8,i)
 790         Dtilde(2,2,i)=b(6,i)
 791         Dtilde(1,1,-i)=b(6,i)
 792         Dtilde(1,2,-i)=-b(8,i)
 793         Dtilde(2,1,-i)=b(8,i)
 794         Dtilde(2,2,-i)=b(6,i)
 795
 796         EEold(1,1,i)= b(10,i)+b(11,i)
 797         eeold(1,1,i)= eenew(1,i)
 798         EEold(2,2,i)=-b(10,i)+b(11,i)
 799         EEold(2,1,i)= b(12,i)-b(13,i)
 800         EEold(1,2,i)= b(12,i)+b(13,i)
 801         EEold(1,1,-i)= b(10,i)+b(11,i)
 802         EEold(2,2,-i)=-b(10,i)+b(11,i)
 803         EEold(2,1,-i)=-b(12,i)+b(13,i)
 804         EEold(1,2,-i)=-b(12,i)-b(13,i)
 805 c        write(iout,*) "TU DOCHODZE"
 806       enddo
 807 c      lprint=.true.
 808       if (lprint) then
 809       do i=1,nloctyp
 810         write (iout,*) 'Type',i
 811         write (iout,*) 'BNEW1(1)'
 812         write (iout,*) (BNEW1(ii,1,i),ii=1,3)
 813         write (iout,*) 'BNEW1(2)'
 814         write (iout,*) BNEW1(1,2,i)
 815         write (iout,*) 'BNEW2(1)'
 816         write (iout,*) (BNEW2(ii,1,i),ii=1,3)
 817         write (iout,*) 'BNEW2(2)'
 818         write (iout,*) BNEW2(1,2,i)
 819         write (iout,*) 'CC'
 820         do j=1,2
 821           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 822         enddo
 823         write(iout,*) 'DD'
 824         do j=1,2
 825           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 826         enddo
 827         write(iout,*) 'EE'
 828         do j=1,2
 829           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
 830         enddo
 831       enddo
 832       endif
 833 #else
 834         B1(1,i)  = b(3,i)
 835         B1(2,i)  = b(5,i)
 836         B1(1,-i) = b(3,i)
 837         B1(2,-i) = -b(5,i)
 838 c        b1(1,i)=0.0d0
 839 c        b1(2,i)=0.0d0
 840         B1tilde(1,i) = b(3,i)
 841         B1tilde(2,i) =-b(5,i)
 842         B1tilde(1,-i) =-b(3,i)
 843         B1tilde(2,-i) =b(5,i)
 844 c        b1tilde(1,i)=0.0d0
 845 c        b1tilde(2,i)=0.0d0
 846         B2(1,i)  = b(2,i)
 847         B2(2,i)  = b(4,i)
 848         B2(1,-i)  =b(2,i)
 849         B2(2,-i)  =-b(4,i)
 850
 851 c        b2(1,i)=0.0d0
 852 c        b2(2,i)=0.0d0
 853         CC(1,1,i)= b(7,i)
 854         CC(2,2,i)=-b(7,i)
 855         CC(2,1,i)= b(9,i)
 856         CC(1,2,i)= b(9,i)
 857         CC(1,1,-i)= b(7,i)
 858         CC(2,2,-i)=-b(7,i)
 859         CC(2,1,-i)=-b(9,i)
 860         CC(1,2,-i)=-b(9,i)
 861 c        CC(1,1,i)=0.0d0
 862 c        CC(2,2,i)=0.0d0
 863 c        CC(2,1,i)=0.0d0
 864 c        CC(1,2,i)=0.0d0
 865         Ctilde(1,1,i)=b(7,i)
 866         Ctilde(1,2,i)=b(9,i)
 867         Ctilde(2,1,i)=-b(9,i)
 868         Ctilde(2,2,i)=b(7,i)
 869         Ctilde(1,1,-i)=b(7,i)
 870         Ctilde(1,2,-i)=-b(9,i)
 871         Ctilde(2,1,-i)=b(9,i)
 872         Ctilde(2,2,-i)=b(7,i)
 873
 874 c        Ctilde(1,1,i)=0.0d0
 875 c        Ctilde(1,2,i)=0.0d0
 876 c        Ctilde(2,1,i)=0.0d0
 877 c        Ctilde(2,2,i)=0.0d0
 878         DD(1,1,i)= b(6,i)
 879         DD(2,2,i)=-b(6,i)
 880         DD(2,1,i)= b(8,i)
 881         DD(1,2,i)= b(8,i)
 882         DD(1,1,-i)= b(6,i)
 883         DD(2,2,-i)=-b(6,i)
 884         DD(2,1,-i)=-b(8,i)
 885         DD(1,2,-i)=-b(8,i)
 886 c        DD(1,1,i)=0.0d0
 887 c        DD(2,2,i)=0.0d0
 888 c        DD(2,1,i)=0.0d0
 889 c        DD(1,2,i)=0.0d0
 890         Dtilde(1,1,i)=b(6,i)
 891         Dtilde(1,2,i)=b(8,i)
 892         Dtilde(2,1,i)=-b(8,i)
 893         Dtilde(2,2,i)=b(6,i)
 894         Dtilde(1,1,-i)=b(6,i)
 895         Dtilde(1,2,-i)=-b(8,i)
 896         Dtilde(2,1,-i)=b(8,i)
 897         Dtilde(2,2,-i)=b(6,i)
 898
 899 c        Dtilde(1,1,i)=0.0d0
 900 c        Dtilde(1,2,i)=0.0d0
 901 c        Dtilde(2,1,i)=0.0d0
 902 c        Dtilde(2,2,i)=0.0d0
 903         EE(1,1,i)= b(10,i)+b(11,i)
 904         EE(2,2,i)=-b(10,i)+b(11,i)
 905         EE(2,1,i)= b(12,i)-b(13,i)
 906         EE(1,2,i)= b(12,i)+b(13,i)
 907         EE(1,1,-i)= b(10,i)+b(11,i)
 908         EE(2,2,-i)=-b(10,i)+b(11,i)
 909         EE(2,1,-i)=-b(12,i)+b(13,i)
 910         EE(1,2,-i)=-b(12,i)-b(13,i)
 911
 912 c        ee(1,1,i)=1.0d0
 913 c        ee(2,2,i)=1.0d0
 914 c        ee(2,1,i)=0.0d0
 915 c        ee(1,2,i)=0.0d0
 916 c        ee(2,1,i)=ee(1,2,i)
 917
 918       enddo
 919       if (lprint) then
 920       do i=1,nloctyp
 921         write (iout,*) 'Type',i
 922         write (iout,*) 'B1'
 923 c        write (iout,'(f10.5)') B1(:,i)
 924         write(iout,*) B1(1,i),B1(2,i)
 925         write (iout,*) 'B2'
 926 c        write (iout,'(f10.5)') B2(:,i)
 927         write(iout,*) B2(1,i),B2(2,i)
 928         write (iout,*) 'CC'
 929         do j=1,2
 930           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 931         enddo
 932         write(iout,*) 'DD'
 933         do j=1,2
 934           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 935         enddo
 936         write(iout,*) 'EE'
 937         do j=1,2
 938           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 939         enddo
 940       enddo
 941       endif
 942 #endif
 943 C
 944 C Read electrostatic-interaction parameters
 945 C
 946       if (lprint) then
 947         write (iout,'(/a)') 'Electrostatic interaction constants:'
 948         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 949      &            'IT','JT','APP','BPP','AEL6','AEL3'
 950       endif
 951       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 952       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 953       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 954       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 955       close (ielep)
 956       do i=1,2
 957         do j=1,2
 958         rri=rpp(i,j)**6
 959         app (i,j)=epp(i,j)*rri*rri
 960         bpp (i,j)=-2.0D0*epp(i,j)*rri
 961         ael6(i,j)=elpp6(i,j)*4.2D0**6
 962         ael3(i,j)=elpp3(i,j)*4.2D0**3
 963         lprint=.true.
 964         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 965      &                    ael6(i,j),ael3(i,j)
 966          lprint=.false.
 967         enddo
 968       enddo
 969 C
 970 C Read side-chain interaction parameters.
 971 C
 972       read (isidep,*) ipot,expon
 973       if (ipot.lt.1 .or. ipot.gt.5) then
 974         write (iout,'(2a)') 'Error while reading SC interaction',
 975      &               'potential file - unknown potential type.'
 976         stop
 977       endif
 978       expon2=expon/2
 979       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 980      & ', exponents are ',expon,2*expon
 981       goto (10,20,30,30,40) ipot
 982 C----------------------- LJ potential ---------------------------------
 983    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 984       if (lprint) then
 985         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 986         write (iout,'(a/)') 'The epsilon array:'
 987         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 988         write (iout,'(/a)') 'One-body parameters:'
 989         write (iout,'(a,4x,a)') 'residue','sigma'
 990         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 991       endif
 992       goto 50
 993 C----------------------- LJK potential --------------------------------
 994    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 995      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 996       if (lprint) then
 997         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 998         write (iout,'(a/)') 'The epsilon array:'
 999         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1000         write (iout,'(/a)') 'One-body parameters:'
1001         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1002         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1003      &        i=1,ntyp)
1004       endif
1005       goto 50
1006 C---------------------- GB or BP potential -----------------------------
1007    30 do i=1,ntyp
1008        read (isidep,*)(eps(i,j),j=i,ntyp)
1009       enddo
1010       read (isidep,*)(sigma0(i),i=1,ntyp)
1011       read (isidep,*)(sigii(i),i=1,ntyp)
1012       read (isidep,*)(chip(i),i=1,ntyp)
1013       read (isidep,*)(alp(i),i=1,ntyp)
1014 C For the GB potential convert sigma'**2 into chi'
1015       if (ipot.eq.4) then
1016         do i=1,ntyp
1017           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1018         enddo
1019       endif
1020       if (lprint) then
1021         write (iout,'(/a/)') 'Parameters of the BP potential:'
1022         write (iout,'(a/)') 'The epsilon array:'
1023         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1024         write (iout,'(/a)') 'One-body parameters:'
1025         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1026      &       '    chip  ','    alph  '
1027         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1028      &                     chip(i),alp(i),i=1,ntyp)
1029       endif
1030       goto 50
1031 C--------------------- GBV potential -----------------------------------
1032    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1033      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1034      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1035       if (lprint) then
1036         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1037         write (iout,'(a/)') 'The epsilon array:'
1038         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1039         write (iout,'(/a)') 'One-body parameters:'
1040         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1041      &      's||/s_|_^2','    chip  ','    alph  '
1042         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1043      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1044       endif
1045    50 continue
1046       close (isidep)
1047 C-----------------------------------------------------------------------
1048 C Calculate the "working" parameters of SC interactions.
1049       do i=2,ntyp
1050         do j=1,i-1
1051           eps(i,j)=eps(j,i)
1052         enddo
1053       enddo
1054       do i=1,ntyp
1055         do j=i,ntyp
1056           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1057           sigma(j,i)=sigma(i,j)
1058           rs0(i,j)=dwa16*sigma(i,j)
1059           rs0(j,i)=rs0(i,j)
1060         enddo
1061       enddo
1062       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1063      & 'Working parameters of the SC interactions:',
1064      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1065      & '  chi1   ','   chi2   '
1066       do i=1,ntyp
1067         do j=i,ntyp
1068           epsij=eps(i,j)
1069           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1070             rrij=sigma(i,j)
1071           else
1072             rrij=rr0(i)+rr0(j)
1073           endif
1074           r0(i,j)=rrij
1075           r0(j,i)=rrij
1076           rrij=rrij**expon
1077           epsij=eps(i,j)
1078           sigeps=dsign(1.0D0,epsij)
1079           epsij=dabs(epsij)
1080           aa(i,j)=epsij*rrij*rrij
1081           bb(i,j)=-sigeps*epsij*rrij
1082           aa(j,i)=aa(i,j)
1083           bb(j,i)=bb(i,j)
1084           if (ipot.gt.2) then
1085             sigt1sq=sigma0(i)**2
1086             sigt2sq=sigma0(j)**2
1087             sigii1=sigii(i)
1088             sigii2=sigii(j)
1089             ratsig1=sigt2sq/sigt1sq
1090             ratsig2=1.0D0/ratsig1
1091             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1092             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1093             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1094           else
1095             rsum_max=sigma(i,j)
1096           endif
1097 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1098             sigmaii(i,j)=rsum_max
1099             sigmaii(j,i)=rsum_max
1100 c         else
1101 c           sigmaii(i,j)=r0(i,j)
1102 c           sigmaii(j,i)=r0(i,j)
1103 c         endif
1104 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1105           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1106             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1107             augm(i,j)=epsij*r_augm**(2*expon)
1108 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1109             augm(j,i)=augm(i,j)
1110           else
1111             augm(i,j)=0.0D0
1112             augm(j,i)=0.0D0
1113           endif
1114           if (lprint) then
1115             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1116      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1117      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1118           endif
1119         enddo
1120       enddo
1121 C
1122 C Define the SC-p interaction constants
1123 C
1124 #ifdef OLDSCP
1125       do i=1,20
1126 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1127 C helix formation)
1128 c       aad(i,1)=0.3D0*4.0D0**12
1129 C Following line for constants currently implemented
1130 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1131         aad(i,1)=1.5D0*4.0D0**12
1132 c       aad(i,1)=0.17D0*5.6D0**12
1133         aad(i,2)=aad(i,1)
1134 C "Soft" SC-p repulsion
1135         bad(i,1)=0.0D0
1136 C Following line for constants currently implemented
1137 c       aad(i,1)=0.3D0*4.0D0**6
1138 C "Hard" SC-p repulsion
1139         bad(i,1)=3.0D0*4.0D0**6
1140 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1141         bad(i,2)=bad(i,1)
1142 c       aad(i,1)=0.0D0
1143 c       aad(i,2)=0.0D0
1144 c       bad(i,1)=1228.8D0
1145 c       bad(i,2)=1228.8D0
1146       enddo
1147 #else
1148 C
1149 C 8/9/01 Read the SC-p interaction constants from file
1150 C
1151       do i=1,ntyp
1152         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1153       enddo
1154       do i=1,ntyp
1155         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1156         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1157         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1158         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1159       enddo
1160
1161       if (lprint) then
1162         write (iout,*) "Parameters of SC-p interactions:"
1163         do i=1,20
1164           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1165      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1166         enddo
1167       endif
1168 #endif
1169 C
1170 C Define the constants of the disulfide bridge
1171 C
1172       ebr=-5.50D0
1173 c
1174 c Old arbitrary potential - commented out.
1175 c
1176 c      dbr= 4.20D0
1177 c      fbr= 3.30D0
1178 c
1179 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1180 c energy surface of diethyl disulfide.
1181 c A. Liwo and U. Kozlowska, 11/24/03
1182 c
1183       D0CM = 3.78d0
1184       AKCM = 15.1d0
1185       AKTH = 11.0d0
1186       AKCT = 12.0d0
1187       V1SS =-1.08d0
1188       V2SS = 7.61d0
1189       V3SS = 13.7d0
1190
1191       if (lprint) then
1192       write (iout,'(/a)') "Disulfide bridge parameters:"
1193       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1194       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1195       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1196       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1197      & ' v3ss:',v3ss
1198       endif
1199       end