source/wham/src-M/parmread.F

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