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