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