source/wham/src-M/parmread.F

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