source/wham/src/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=-ntyp1,-1
 303         ithetyp(i)=-ithetyp(-i)
 304       enddo
 305 c      write (iout,*) "tu dochodze"
 306       do iblock=1,2
 307       do i=-maxthetyp,maxthetyp
 308         do j=-maxthetyp,maxthetyp
 309           do k=-maxthetyp,maxthetyp
 310             aa0thet(i,j,k,iblock)=0.0d0
 311             do l=1,ntheterm
 312               aathet(l,i,j,k,iblock)=0.0d0
 313             enddo
 314             do l=1,ntheterm2
 315               do m=1,nsingle
 316                 bbthet(m,l,i,j,k,iblock)=0.0d0
 317                 ccthet(m,l,i,j,k,iblock)=0.0d0
 318                 ddthet(m,l,i,j,k,iblock)=0.0d0
 319                 eethet(m,l,i,j,k,iblock)=0.0d0
 320               enddo
 321             enddo
 322             do l=1,ntheterm3
 323               do m=1,ndouble
 324                 do mm=1,ndouble
 325                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 326                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 327                 enddo
 328               enddo
 329             enddo
 330           enddo
 331         enddo
 332       enddo
 333       enddo
 334
 335       do iblock=1,2
 336       do i=0,nthetyp
 337         do j=-nthetyp,nthetyp
 338           do k=-nthetyp,nthetyp
 339             read (ithep,'(6a)') res1
 340             read (ithep,*) aa0thet(i,j,k,iblock)
 341             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 342             read (ithep,*)
 343      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 344      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 345      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 346      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 347      &        ,ll=1,ntheterm2)
 348             read (ithep,*)
 349      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 350      &      ffthet(lll,llll,ll,i,j,k,iblock),
 351      &         ggthet(llll,lll,ll,i,j,k,iblock)
 352      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 353      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 354           enddo
 355         enddo
 356       enddo
 357 C
 358 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 359 C coefficients of theta-and-gamma-dependent terms are zero.
 360 C
 361       do i=1,nthetyp
 362         do j=1,nthetyp
 363           do l=1,ntheterm
 364             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 365             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 366           enddo
 367           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 368           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 369         enddo
 370         do l=1,ntheterm
 371           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 372         enddo
 373         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 374       enddo
 375       enddo
 376 C Substitution for D aminoacids from symmetry.
 377       do iblock=1,2
 378       do i=-nthetyp,0
 379         do j=-nthetyp,nthetyp
 380           do k=-nthetyp,nthetyp
 381            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 382            do l=1,ntheterm
 383            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 384            enddo
 385            do ll=1,ntheterm2
 386             do lll=1,nsingle
 387             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 388             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 389             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 390             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 391             enddo
 392           enddo
 393           do ll=1,ntheterm3
 394            do lll=2,ndouble
 395             do llll=1,lll-1
 396             ffthet(llll,lll,ll,i,j,k,iblock)=
 397      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 398             ffthet(lll,llll,ll,i,j,k,iblock)=
 399      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 400             ggthet(llll,lll,ll,i,j,k,iblock)=
 401      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 402             ggthet(lll,llll,ll,i,j,k,iblock)=
 403      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 404             enddo !ll
 405            enddo  !lll
 406           enddo   !llll
 407          enddo    !k
 408         enddo     !j
 409        enddo      !i
 410       enddo       !iblock
 411 C
 412 C Control printout of the coefficients of virtual-bond-angle potentials
 413 C
 414       if (lprint) then
 415         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 416         do i=1,nthetyp+1
 417           do j=1,nthetyp+1
 418             do k=1,nthetyp+1
 419               write (iout,'(//4a)')
 420      &         'Type ',onelett(i),onelett(j),onelett(k)
 421               write (iout,'(//a,10x,a)') " l","a[l]"
 422               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 423               write (iout,'(i2,1pe15.5)')
 424      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 425             do l=1,ntheterm2
 426               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 427      &          "b",l,"c",l,"d",l,"e",l
 428               do m=1,nsingle
 429                 write (iout,'(i2,4(1pe15.5))') m,
 430      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 431      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 432               enddo
 433             enddo
 434             do l=1,ntheterm3
 435               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 436      &          "f+",l,"f-",l,"g+",l,"g-",l
 437               do m=2,ndouble
 438                 do n=1,m-1
 439                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 440      &              ffthet(n,m,l,i,j,k,iblock),
 441      &              ffthet(m,n,l,i,j,k,iblock),
 442      &              ggthet(n,m,l,i,j,k,iblock),
 443      &              ggthet(m,n,l,i,j,k,iblock)
 444                 enddo
 445               enddo
 446             enddo
 447           enddo
 448         enddo
 449       enddo
 450       call flush(iout)
 451       endif
 452 #endif
 453
 454 #ifdef CRYST_SC
 455 C
 456 C Read the parameters of the probability distribution/energy expression
 457 C of the side chains.
 458 C
 459       do i=1,ntyp
 460         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 461         if (i.eq.10) then
 462           dsc_inv(i)=0.0D0
 463         else
 464           dsc_inv(i)=1.0D0/dsc(i)
 465         endif
 466         if (i.ne.10) then
 467         do j=1,nlob(i)
 468           do k=1,3
 469             do l=1,3
 470               blower(l,k,j)=0.0D0
 471             enddo
 472           enddo
 473         enddo
 474         bsc(1,i)=0.0D0
 475         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 476          censc(1,1,-i)=censc(1,1,i)
 477          censc(2,1,-i)=censc(2,1,i)
 478          censc(3,1,-i)=-censc(3,1,i)
 479         do j=2,nlob(i)
 480           read (irotam,*) bsc(j,i)
 481           read (irotam,*) (censc(k,j,i),k=1,3),
 482      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 483           censc(1,j,-i)=censc(1,j,i)
 484           censc(2,j,-i)=censc(2,j,i)
 485           censc(3,j,-i)=-censc(3,j,i)
 486         enddo
 487         do j=1,nlob(i)
 488           do k=1,3
 489             do l=1,k
 490               akl=0.0D0
 491               do m=1,3
 492                 akl=akl+blower(k,m,j)*blower(l,m,j)
 493               enddo
 494               gaussc(k,l,j,i)=akl
 495               gaussc(l,k,j,i)=akl
 496                if (((k.eq.3).and.(l.ne.3))
 497      &        .or.((l.eq.3).and.(k.ne.3))) then
 498                  gaussc(k,l,j,-i)=-akl
 499                  gaussc(l,k,j,-i)=-akl
 500                else
 501                  gaussc(k,l,j,-i)=akl
 502                  gaussc(l,k,j,-i)=akl
 503                endif
 504             enddo
 505           enddo
 506         enddo
 507         endif
 508       enddo
 509       close (irotam)
 510       if (lprint) then
 511         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 512         do i=1,ntyp
 513           nlobi=nlob(i)
 514           if (nlobi.gt.0) then
 515           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 516      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 517 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 518 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 519 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 520 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 521            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 522      &                             'log h',(bsc(j,i),j=1,nlobi)
 523            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 524      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 525 c          write (iout,'(a)')
 526 c         do j=1,nlobi
 527 c           ind=0
 528 c           do k=1,3
 529 c             do l=1,k
 530 c              ind=ind+1
 531 c              blower(k,l,j)=gaussc(ind,j,i)
 532 c             enddo
 533 c           enddo
 534 c         enddo
 535           do k=1,3
 536             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 537      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 538           enddo
 539           endif
 540         enddo
 541       endif
 542 #else
 543 C
 544 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 545 C added by Urszula Kozlowska 07/11/2007
 546 C
 547       do i=1,ntyp
 548         read (irotam,*)
 549        if (i.eq.10) then
 550          read (irotam,*)
 551        else
 552          do j=1,65
 553            read(irotam,*) sc_parmin(j,i)
 554          enddo
 555        endif
 556       enddo
 557 #endif
 558       close(irotam)
 559 #ifdef CRYST_TOR
 560 C
 561 C Read torsional parameters in old format
 562 C
 563       read (itorp,*) ntortyp,nterm_old
 564       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 565       read (itorp,*) (itortyp(i),i=1,ntyp)
 566       do i=1,ntortyp
 567         do j=1,ntortyp
 568           read (itorp,'(a)')
 569           do k=1,nterm_old
 570             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 571           enddo
 572         enddo
 573       enddo
 574       close (itorp)
 575       if (lprint) then
 576         write (iout,'(/a/)') 'Torsional constants:'
 577         do i=1,ntortyp
 578           do j=1,ntortyp
 579             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 580             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 581           enddo
 582         enddo
 583       endif
 584
 585
 586 #else
 587 C
 588 C Read torsional parameters
 589 C
 590       read (itorp,*) ntortyp
 591       read (itorp,*) (itortyp(i),i=1,ntyp)
 592       do iblock=1,2
 593       do i=-ntyp,-1
 594        itortyp(i)=-itortyp(-i)
 595       enddo
 596 c      write (iout,*) 'ntortyp',ntortyp
 597       do i=0,ntortyp-1
 598         do j=-ntortyp+1,ntortyp-1
 599           read (itorp,*) nterm(i,j,iblock),
 600      &          nlor(i,j,iblock)
 601           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 602           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 603           v0ij=0.0d0
 604           si=-1.0d0
 605           do k=1,nterm(i,j,iblock)
 606             read (itorp,*) kk,v1(k,i,j,iblock),v2(k,i,j,iblock)
 607             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 608             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 609             v0ij=v0ij+si*v1(k,i,j,iblock)
 610             si=-si
 611           enddo
 612           do k=1,nlor(i,j,iblock)
 613             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 614             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 615           enddo
 616           v0(i,j,iblock)=v0ij
 617           v0(-i,-j,iblock)=v0ij
 618         enddo
 619       enddo
 620       enddo
 621       close (itorp)
 622       if (lprint) then
 623         write (iout,'(/a/)') 'Torsional constants:'
 624         do i=1,ntortyp
 625           do j=1,ntortyp
 626             write (iout,*) 'ityp',i,' jtyp',j
 627             write (iout,*) 'Fourier constants'
 628             do k=1,nterm(i,j,iblock)
 629               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 630      &        v2(k,i,j,iblock)
 631             enddo
 632             write (iout,*) 'Lorenz constants'
 633             do k=1,nlor(i,j,iblock)
 634               write (iout,'(3(1pe15.5))')
 635      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 636             enddo
 637           enddo
 638         enddo
 639       endif
 640 C
 641 C 6/23/01 Read parameters for double torsionals
 642 C
 643       do iblock=1,2
 644       do i=0,ntortyp-1
 645         do j=-ntortyp+1,ntortyp-1
 646           do k=-ntortyp+1,ntortyp-1
 647             read (itordp,'(3a1)') t1,t2,t3
 648             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 649      &        .or. t3.ne.toronelet(k)) then
 650               write (iout,*) "Error in double torsional parameter file",
 651      &         i,j,k,t1,t2,t3
 652                stop "Error in double torsional parameter file"
 653             endif
 654          read (itordp,*) ntermd_1(i,j,k,iblock),
 655      &         ntermd_2(i,j,k,iblock)
 656             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 657             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 658             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
 659      &         ntermd_1(i,j,k,iblock))
 660             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
 661      &         ntermd_1(i,j,k,iblock))
 662             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
 663      &         ntermd_1(i,j,k,iblock))
 664             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
 665      &         ntermd_1(i,j,k,iblock))
 666 C Martix of D parameters for one dimesional foureir series
 667             do l=1,ntermd_1(i,j,k,iblock)
 668              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 669              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 670              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 671              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 672 c            write(iout,*) "whcodze" ,
 673 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 674             enddo
 675             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
 676      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 677      &         v2s(m,l,i,j,k,iblock),
 678      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 679 C Martix of D parameters for two dimesional fourier series
 680             do l=1,ntermd_2(i,j,k,iblock)
 681              do m=1,l-1
 682              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 683              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 684              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 685              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 686           enddo!m
 687         enddo!l
 688       enddo!k
 689       enddo!j
 690       enddo!i
 691       enddo!iblock
 692       if (lprint) then
 693       write (iout,*)
 694       write (iout,*) 'Constants for double torsionals'
 695       do iblock=1,2
 696       do i=0,ntortyp-1
 697         do j=-ntortyp+1,ntortyp-1
 698           do k=-ntortyp+1,ntortyp-1
 699             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 700      &        ' nsingle',ntermd_1(i,j,k,iblock),
 701      &        ' ndouble',ntermd_2(i,j,k,iblock)
 702             write (iout,*)
 703             write (iout,*) 'Single angles:'
 704             do l=1,ntermd_1(i,j,k,iblock)
 705               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 706      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 707      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
 708             enddo
 709             write (iout,*)
 710             write (iout,*) 'Pairs of angles:'
 711             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 712             do l=1,ntermd_2(i,j,k,iblock)
 713               write (iout,'(i5,20f10.5)')
 714      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 715             enddo
 716             write (iout,*)
 717             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 718             do l=1,ntermd_2(i,j,k,iblock)
 719               write (iout,'(i5,20f10.5)')
 720      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 721             enddo
 722             write (iout,*)
 723           enddo
 724         enddo
 725       enddo
 726       enddo
 727       endif
 728 #endif
 729 C Read of Side-chain backbone correlation parameters
 730 C Modified 11 May 2012 by Adasko
 731 CCC
 732 C
 733       read (isccor,*) nsccortyp
 734       read (isccor,*) (isccortyp(i),i=1,ntyp)
 735       do i=-ntyp,-1
 736         isccortyp(i)=-isccortyp(-i)
 737        enddo
 738        iscprol=isccortyp(20)
 739 c      write (iout,*) 'ntortyp',ntortyp
 740       maxinter=3
 741 cc maxinter is maximum interaction sites
 742       do l=1,maxinter
 743       do i=1,nsccortyp
 744         do j=1,nsccortyp
 745           read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
 746           v0ijsccor=0.0d0
 747           v0ijsccor1=0.0d0
 748           v0ijsccor2=0.0d0
 749           v0ijsccor3=0.0d0
 750           si=-1.0d0
 751           nterm_sccor(-i,j)=nterm_sccor(i,j)
 752           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 753           nterm_sccor(i,-j)=nterm_sccor(i,j)
 754           do k=1,nterm_sccor(i,j)
 755             read (isccor,*) kk,v1sccor(k,l,i,j)
 756      &    ,v2sccor(k,l,i,j)
 757              if (j.eq.iscprol) then
 758              if (i.eq.isccortyp(10)) then
 759              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 760              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 761              else
 762               v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 763      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 764              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 765      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 766              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 767              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 768              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 769              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 770              endif
 771             else
 772              if (i.eq.isccortyp(10)) then
 773              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 774              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 775              else
 776                if (j.eq.isccortyp(10)) then
 777              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 778              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 779                else
 780              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 781              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 782              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 783              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 784              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 785              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 786                 endif
 787                endif
 788              endif
 789             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 790             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 791             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 792             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 793             si=-si
 794           enddo
 795           do k=1,nlor_sccor(i,j)
 796             read (isccor,*) kk,vlor1sccor(k,i,j),
 797      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 798             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 799      &(1+vlor3sccor(k,i,j)**2)
 800           enddo
 801           v0sccor(l,i,j)=v0ijsccor
 802           v0sccor(l,-i,j)=v0ijsccor1
 803           v0sccor(l,i,-j)=v0ijsccor2
 804           v0sccor(l,-i,-j)=v0ijsccor3
 805         enddo
 806       enddo
 807       enddo
 808       close (isccor)
 809
 810       if (lprint) then
 811         write (iout,'(/a/)') 'Torsional constants:'
 812         do i=1,nsccortyp
 813           do j=1,nsccortyp
 814             write (iout,*) 'ityp',i,' jtyp',j
 815             write (iout,*) 'Fourier constants'
 816             do k=1,nterm_sccor(i,j)
 817            write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 818             enddo
 819             write (iout,*) 'Lorenz constants'
 820             do k=1,nlor_sccor(i,j)
 821               write (iout,'(3(1pe15.5))')
 822      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 823             enddo
 824           enddo
 825         enddo
 826       endif
 827
 828 C
 829 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 830 C         interaction energy of the Gly, Ala, and Pro prototypes.
 831 C
 832       read (ifourier,*) nloctyp
 833       do i=0,nloctyp-1
 834         read (ifourier,*)
 835         read (ifourier,*) (b(ii),ii=1,13)
 836         if (lprint) then
 837         write (iout,*) 'Type',i
 838         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 839         endif
 840         B1(1,i)  = b(3)
 841         B1(2,i)  = b(5)
 842         B1(1,-i) = b(3)
 843         B1(2,-i) = -b(5)
 844 c        b1(1,i)=0.0d0
 845 c        b1(2,i)=0.0d0
 846         B1tilde(1,i) = b(3)
 847         B1tilde(2,i) =-b(5)
 848         B1tilde(1,-i) =-b(3)
 849         B1tilde(2,-i) =b(5)
 850 c        b1tilde(1,i)=0.0d0
 851 c        b1tilde(2,i)=0.0d0
 852         B2(1,i)  = b(2)
 853         B2(2,i)  = b(4)
 854         B2(1,-i)  =b(2)
 855         B2(2,-i)  =-b(4)
 856
 857 c        b2(1,i)=0.0d0
 858 c        b2(2,i)=0.0d0
 859         CC(1,1,i)= b(7)
 860         CC(2,2,i)=-b(7)
 861         CC(2,1,i)= b(9)
 862         CC(1,2,i)= b(9)
 863         CC(1,1,-i)= b(7)
 864         CC(2,2,-i)=-b(7)
 865         CC(2,1,-i)=-b(9)
 866         CC(1,2,-i)=-b(9)
 867 c        CC(1,1,i)=0.0d0
 868 c        CC(2,2,i)=0.0d0
 869 c        CC(2,1,i)=0.0d0
 870 c        CC(1,2,i)=0.0d0
 871         Ctilde(1,1,i)=b(7)
 872         Ctilde(1,2,i)=b(9)
 873         Ctilde(2,1,i)=-b(9)
 874         Ctilde(2,2,i)=b(7)
 875         Ctilde(1,1,-i)=b(7)
 876         Ctilde(1,2,-i)=-b(9)
 877         Ctilde(2,1,-i)=b(9)
 878         Ctilde(2,2,-i)=b(7)
 879
 880 c        Ctilde(1,1,i)=0.0d0
 881 c        Ctilde(1,2,i)=0.0d0
 882 c        Ctilde(2,1,i)=0.0d0
 883 c        Ctilde(2,2,i)=0.0d0
 884         DD(1,1,i)= b(6)
 885         DD(2,2,i)=-b(6)
 886         DD(2,1,i)= b(8)
 887         DD(1,2,i)= b(8)
 888         DD(1,1,-i)= b(6)
 889         DD(2,2,-i)=-b(6)
 890         DD(2,1,-i)=-b(8)
 891         DD(1,2,-i)=-b(8)
 892 c        DD(1,1,i)=0.0d0
 893 c        DD(2,2,i)=0.0d0
 894 c        DD(2,1,i)=0.0d0
 895 c        DD(1,2,i)=0.0d0
 896         Dtilde(1,1,i)=b(6)
 897         Dtilde(1,2,i)=b(8)
 898         Dtilde(2,1,i)=-b(8)
 899         Dtilde(2,2,i)=b(6)
 900         Dtilde(1,1,-i)=b(6)
 901         Dtilde(1,2,-i)=-b(8)
 902         Dtilde(2,1,-i)=b(8)
 903         Dtilde(2,2,-i)=b(6)
 904
 905 c        Dtilde(1,1,i)=0.0d0
 906 c        Dtilde(1,2,i)=0.0d0
 907 c        Dtilde(2,1,i)=0.0d0
 908 c        Dtilde(2,2,i)=0.0d0
 909         EE(1,1,i)= b(10)+b(11)
 910         EE(2,2,i)=-b(10)+b(11)
 911         EE(2,1,i)= b(12)-b(13)
 912         EE(1,2,i)= b(12)+b(13)
 913         EE(1,1,-i)= b(10)+b(11)
 914         EE(2,2,-i)=-b(10)+b(11)
 915         EE(2,1,-i)=-b(12)+b(13)
 916         EE(1,2,-i)=-b(12)-b(13)
 917
 918 c        ee(1,1,i)=1.0d0
 919 c        ee(2,2,i)=1.0d0
 920 c        ee(2,1,i)=0.0d0
 921 c        ee(1,2,i)=0.0d0
 922 c        ee(2,1,i)=ee(1,2,i)
 923       enddo
 924       if (lprint) then
 925       do i=1,nloctyp
 926         write (iout,*) 'Type',i
 927         write (iout,*) 'B1'
 928 c        write (iout,'(f10.5)') B1(:,i)
 929         write(iout,*) B1(1,i),B1(2,i)
 930         write (iout,*) 'B2'
 931 c        write (iout,'(f10.5)') B2(:,i)
 932         write(iout,*) B2(1,i),B2(2,i)
 933         write (iout,*) 'CC'
 934         do j=1,2
 935           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 936         enddo
 937         write(iout,*) 'DD'
 938         do j=1,2
 939           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 940         enddo
 941         write(iout,*) 'EE'
 942         do j=1,2
 943           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 944         enddo
 945       enddo
 946       endif
 947 C
 948 C Read electrostatic-interaction parameters
 949 C
 950       if (lprint) then
 951         write (iout,'(/a)') 'Electrostatic interaction constants:'
 952         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 953      &            'IT','JT','APP','BPP','AEL6','AEL3'
 954       endif
 955       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 956       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 957       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 958       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 959       close (ielep)
 960       do i=1,2
 961         do j=1,2
 962         rri=rpp(i,j)**6
 963         app (i,j)=epp(i,j)*rri*rri
 964         bpp (i,j)=-2.0D0*epp(i,j)*rri
 965         ael6(i,j)=elpp6(i,j)*4.2D0**6
 966         ael3(i,j)=elpp3(i,j)*4.2D0**3
 967         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 968      &                    ael6(i,j),ael3(i,j)
 969         enddo
 970       enddo
 971 C
 972 C Read side-chain interaction parameters.
 973 C
 974       read (isidep,*) ipot,expon
 975       if (ipot.lt.1 .or. ipot.gt.5) then
 976         write (iout,'(2a)') 'Error while reading SC interaction',
 977      &               'potential file - unknown potential type.'
 978         stop
 979       endif
 980       expon2=expon/2
 981       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 982      & ', exponents are ',expon,2*expon
 983       goto (10,20,30,30,40) ipot
 984 C----------------------- LJ potential ---------------------------------
 985    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 986       if (lprint) then
 987         write (iout,'(/a/)') 'Parameters of the LJ 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,a)') 'residue','sigma'
 992         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 993       endif
 994       goto 50
 995 C----------------------- LJK potential --------------------------------
 996    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 997      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 998       if (lprint) then
 999         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1000         write (iout,'(a/)') 'The epsilon array:'
1001         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1002         write (iout,'(/a)') 'One-body parameters:'
1003         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1004         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1005      &        i=1,ntyp)
1006       endif
1007       goto 50
1008 C---------------------- GB or BP potential -----------------------------
1009    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1010      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
1011      &  (alp(i),i=1,ntyp)
1012 C For the GB potential convert sigma'**2 into chi'
1013       if (ipot.eq.4) then
1014         do i=1,ntyp
1015           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1016         enddo
1017       endif
1018       if (lprint) then
1019         write (iout,'(/a/)') 'Parameters of the BP potential:'
1020         write (iout,'(a/)') 'The epsilon array:'
1021         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1022         write (iout,'(/a)') 'One-body parameters:'
1023         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1024      &       '    chip  ','    alph  '
1025         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1026      &                     chip(i),alp(i),i=1,ntyp)
1027       endif
1028       goto 50
1029 C--------------------- GBV potential -----------------------------------
1030    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1031      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1032      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1033       if (lprint) then
1034         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1035         write (iout,'(a/)') 'The epsilon array:'
1036         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1037         write (iout,'(/a)') 'One-body parameters:'
1038         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1039      &      's||/s_|_^2','    chip  ','    alph  '
1040         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1041      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1042       endif
1043    50 continue
1044       close (isidep)
1045 C-----------------------------------------------------------------------
1046 C Calculate the "working" parameters of SC interactions.
1047       do i=2,ntyp
1048         do j=1,i-1
1049           eps(i,j)=eps(j,i)
1050         enddo
1051       enddo
1052       do i=1,ntyp
1053         do j=i,ntyp
1054           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1055           sigma(j,i)=sigma(i,j)
1056           rs0(i,j)=dwa16*sigma(i,j)
1057           rs0(j,i)=rs0(i,j)
1058         enddo
1059       enddo
1060       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1061      & 'Working parameters of the SC interactions:',
1062      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1063      & '  chi1   ','   chi2   '
1064       do i=1,ntyp
1065         do j=i,ntyp
1066           epsij=eps(i,j)
1067           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1068             rrij=sigma(i,j)
1069           else
1070             rrij=rr0(i)+rr0(j)
1071           endif
1072           r0(i,j)=rrij
1073           r0(j,i)=rrij
1074           rrij=rrij**expon
1075           epsij=eps(i,j)
1076           sigeps=dsign(1.0D0,epsij)
1077           epsij=dabs(epsij)
1078           aa(i,j)=epsij*rrij*rrij
1079           bb(i,j)=-sigeps*epsij*rrij
1080           aa(j,i)=aa(i,j)
1081           bb(j,i)=bb(i,j)
1082           if (ipot.gt.2) then
1083             sigt1sq=sigma0(i)**2
1084             sigt2sq=sigma0(j)**2
1085             sigii1=sigii(i)
1086             sigii2=sigii(j)
1087             ratsig1=sigt2sq/sigt1sq
1088             ratsig2=1.0D0/ratsig1
1089             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1090             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1091             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1092           else
1093             rsum_max=sigma(i,j)
1094           endif
1095 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1096             sigmaii(i,j)=rsum_max
1097             sigmaii(j,i)=rsum_max
1098 c         else
1099 c           sigmaii(i,j)=r0(i,j)
1100 c           sigmaii(j,i)=r0(i,j)
1101 c         endif
1102 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1103           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1104             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1105             augm(i,j)=epsij*r_augm**(2*expon)
1106 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1107             augm(j,i)=augm(i,j)
1108           else
1109             augm(i,j)=0.0D0
1110             augm(j,i)=0.0D0
1111           endif
1112           if (lprint) then
1113             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1114      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1115      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1116           endif
1117         enddo
1118       enddo
1119 C
1120 C Define the SC-p interaction constants
1121 C
1122 #ifdef OLDSCP
1123       do i=1,20
1124 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1125 C helix formation)
1126 c       aad(i,1)=0.3D0*4.0D0**12
1127 C Following line for constants currently implemented
1128 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1129         aad(i,1)=1.5D0*4.0D0**12
1130 c       aad(i,1)=0.17D0*5.6D0**12
1131         aad(i,2)=aad(i,1)
1132 C "Soft" SC-p repulsion
1133         bad(i,1)=0.0D0
1134 C Following line for constants currently implemented
1135 c       aad(i,1)=0.3D0*4.0D0**6
1136 C "Hard" SC-p repulsion
1137         bad(i,1)=3.0D0*4.0D0**6
1138 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1139         bad(i,2)=bad(i,1)
1140 c       aad(i,1)=0.0D0
1141 c       aad(i,2)=0.0D0
1142 c       bad(i,1)=1228.8D0
1143 c       bad(i,2)=1228.8D0
1144       enddo
1145 #else
1146 C
1147 C 8/9/01 Read the SC-p interaction constants from file
1148 C
1149       do i=1,ntyp
1150         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1151       enddo
1152       do i=1,ntyp
1153         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1154         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1155         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1156         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1157       enddo
1158
1159       if (lprint) then
1160         write (iout,*) "Parameters of SC-p interactions:"
1161         do i=1,20
1162           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1163      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1164         enddo
1165       endif
1166 #endif
1167 C
1168 C Define the constants of the disulfide bridge
1169 C
1170       ebr=-5.50D0
1171 c
1172 c Old arbitrary potential - commented out.
1173 c
1174 c      dbr= 4.20D0
1175 c      fbr= 3.30D0
1176 c
1177 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1178 c energy surface of diethyl disulfide.
1179 c A. Liwo and U. Kozlowska, 11/24/03
1180 c
1181       D0CM = 3.78d0
1182       AKCM = 15.1d0
1183       AKTH = 11.0d0
1184       AKCT = 12.0d0
1185       V1SS =-1.08d0
1186       V2SS = 7.61d0
1187       V3SS = 13.7d0
1188
1189       if (lprint) then
1190       write (iout,'(/a)') "Disulfide bridge parameters:"
1191       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1192       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1193       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1194       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1195      & ' v3ss:',v3ss
1196       endif
1197       return
1198       end