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