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