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