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