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       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
  55       write (iout,*) "iparm",iparm," myparm",myparm
  56 c If reading not own parameters, skip assignment
  57
  58       if (iparm.eq.myparm .or. .not.separate_parset) then
  59
  60 c
  61 c Setup weights for UNRES
  62 c
  63       wsc=ww(1)
  64       wscp=ww(2)
  65       welec=ww(3)
  66       wcorr=ww(4)
  67       wcorr5=ww(5)
  68       wcorr6=ww(6)
  69       wel_loc=ww(7)
  70       wturn3=ww(8)
  71       wturn4=ww(9)
  72       wturn6=ww(10)
  73       wang=ww(11)
  74       wscloc=ww(12)
  75       wtor=ww(13)
  76       wtor_d=ww(14)
  77       wvdwpp=ww(16)
  78       wbond=ww(18)
  79       wsccor=ww(19)
  80
  81       endif
  82
  83       call card_concat(controlcard,.false.)
  84
  85 c Return if not own parameters
  86
  87       if (iparm.ne.myparm .and. separate_parset) return
  88
  89       call reads(controlcard,"BONDPAR",bondname_t,bondname)
  90       open (ibond,file=bondname_t,status='old')
  91       rewind(ibond)
  92       call reads(controlcard,"THETPAR",thetname_t,thetname)
  93       open (ithep,file=thetname_t,status='old')
  94       rewind(ithep)
  95       call reads(controlcard,"ROTPAR",rotname_t,rotname)
  96       open (irotam,file=rotname_t,status='old')
  97       rewind(irotam)
  98       call reads(controlcard,"TORPAR",torname_t,torname)
  99       open (itorp,file=torname_t,status='old')
 100       rewind(itorp)
 101       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 102       open (itordp,file=tordname_t,status='old')
 103       rewind(itordp)
 104       call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
 105       open (isccor,file=sccorname_t,status='old')
 106       rewind(isccor)
 107       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 108       open (ifourier,file=fouriername_t,status='old')
 109       rewind(ifourier)
 110       call reads(controlcard,"ELEPAR",elename_t,elename)
 111       open (ielep,file=elename_t,status='old')
 112       rewind(ielep)
 113       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 114       open (isidep,file=sidename_t,status='old')
 115       rewind(isidep)
 116       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 117       open (iscpp,file=scpname_t,status='old')
 118       rewind(iscpp)
 119       write (iout,*) "Parameter set:",iparm
 120       write (iout,*) "Energy-term weights:"
 121       do i=1,n_ene
 122         write (iout,'(a16,f10.5)') wname(i),ww(i)
 123       enddo
 124       write (iout,*) "Sidechain potential file        : ",
 125      &  sidename_t(:ilen(sidename_t))
 126 #ifndef OLDSCP
 127       write (iout,*) "SCp potential file              : ",
 128      &  scpname_t(:ilen(scpname_t))
 129 #endif
 130       write (iout,*) "Electrostatic potential file    : ",
 131      &  elename_t(:ilen(elename_t))
 132       write (iout,*) "Cumulant coefficient file       : ",
 133      &  fouriername_t(:ilen(fouriername_t))
 134       write (iout,*) "Torsional parameter file        : ",
 135      &  torname_t(:ilen(torname_t))
 136       write (iout,*) "Double torsional parameter file : ",
 137      &  tordname_t(:ilen(tordname_t))
 138       write (iout,*) "Backbone-rotamer parameter file : ",
 139      &  sccorname(:ilen(sccorname))
 140       write (iout,*) "Bond & inertia constant file    : ",
 141      &  bondname_t(:ilen(bondname_t))
 142       write (iout,*) "Bending parameter file          : ",
 143      &  thetname_t(:ilen(thetname_t))
 144       write (iout,*) "Rotamer parameter file          : ",
 145      &  rotname_t(:ilen(rotname_t))
 146
 147 c
 148 c Read the virtual-bond parameters, masses, and moments of inertia
 149 c and Stokes' radii of the peptide group and side chains
 150 c
 151 #ifdef CRYST_BOND
 152       read (ibond,*) vbldp0,akp
 153       do i=1,ntyp
 154         nbondterm(i)=1
 155         read (ibond,*) vbldsc0(1,i),aksc(1,i)
 156         dsc(i) = vbldsc0(1,i)
 157         if (i.eq.10) then
 158           dsc_inv(i)=0.0D0
 159         else
 160           dsc_inv(i)=1.0D0/dsc(i)
 161         endif
 162       enddo
 163 #else
 164       read (ibond,*) ijunk,vbldp0,akp,rjunk
 165       do i=1,ntyp
 166         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
 167      &   j=1,nbondterm(i))
 168         dsc(i) = vbldsc0(1,i)
 169         if (i.eq.10) then
 170           dsc_inv(i)=0.0D0
 171         else
 172           dsc_inv(i)=1.0D0/dsc(i)
 173         endif
 174       enddo
 175 #endif
 176       if (lprint) then
 177         write(iout,'(/a/)')"Force constants virtual bonds:"
 178         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
 179      &   'inertia','Pstok'
 180         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 181         do i=1,ntyp
 182           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
 183      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
 184           do j=2,nbondterm(i)
 185             write (iout,'(13x,3f10.5)')
 186      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
 187           enddo
 188         enddo
 189       endif
 190 #ifdef CRYST_THETA
 191 C
 192 C Read the parameters of the probability distribution/energy expression
 193 C of the virtual-bond valence angles theta
 194 C
 195       do i=1,ntyp
 196         read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 197         read (ithep,*) (polthet(j,i),j=0,3)
 198         read (ithep,*) (gthet(j,i),j=1,3)
 199         read (ithep,*) theta0(i),sig0(i),sigc0(i)
 200         sigc0(i)=sigc0(i)**2
 201       enddo
 202       close (ithep)
 203       if (lprint) then
 204 c       write (iout,'(a)')
 205 c    &    'Parameters of the virtual-bond valence angles:'
 206 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 207 c    & '    ATHETA0   ','         A1   ','        A2    ',
 208 c    & '        B1    ','         B2   '
 209 c       do i=1,ntyp
 210 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 211 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 212 c       enddo
 213 c       write (iout,'(/a/9x,5a/79(1h-))')
 214 c    & 'Parameters of the expression for sigma(theta_c):',
 215 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 216 c    & '     ALPH3    ','    SIGMA0C   '
 217 c       do i=1,ntyp
 218 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 219 c    &      (polthet(j,i),j=0,3),sigc0(i)
 220 c       enddo
 221 c       write (iout,'(/a/9x,5a/79(1h-))')
 222 c    & 'Parameters of the second gaussian:',
 223 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 224 c    & '        G2    ','         G3   '
 225 c       do i=1,ntyp
 226 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 227 c    &       sig0(i),(gthet(j,i),j=1,3)
 228 c       enddo
 229         write (iout,'(a)')
 230      &    'Parameters of the virtual-bond valence angles:'
 231         write (iout,'(/a/9x,5a/79(1h-))')
 232      & 'Coefficients of expansion',
 233      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 234      & '   b1*10^1    ','    b2*10^1   '
 235         do i=1,ntyp
 236           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 237      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
 238         enddo
 239         write (iout,'(/a/9x,5a/79(1h-))')
 240      & 'Parameters of the expression for sigma(theta_c):',
 241      & ' alpha0       ','  alph1       ',' alph2        ',
 242      & ' alhp3        ','   sigma0c    '
 243         do i=1,ntyp
 244           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 245      &      (polthet(j,i),j=0,3),sigc0(i)
 246         enddo
 247         write (iout,'(/a/9x,5a/79(1h-))')
 248      & 'Parameters of the second gaussian:',
 249      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 250      & '        G2    ','   G3*10^1    '
 251         do i=1,ntyp
 252           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 253      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 254         enddo
 255       endif
 256 #else
 257 C
 258 C Read the parameters of Utheta determined from ab initio surfaces
 259 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 260 C
 261       read (ithep,*) nthetyp,ntheterm,ntheterm2,
 262      &  ntheterm3,nsingle,ndouble
 263       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 264       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 265       do i=1,maxthetyp
 266         do j=1,maxthetyp
 267           do k=1,maxthetyp
 268             aa0thet(i,j,k)=0.0d0
 269             do l=1,ntheterm
 270               aathet(l,i,j,k)=0.0d0
 271             enddo
 272             do l=1,ntheterm2
 273               do m=1,nsingle
 274                 bbthet(m,l,i,j,k)=0.0d0
 275                 ccthet(m,l,i,j,k)=0.0d0
 276                 ddthet(m,l,i,j,k)=0.0d0
 277                 eethet(m,l,i,j,k)=0.0d0
 278               enddo
 279             enddo
 280             do l=1,ntheterm3
 281               do m=1,ndouble
 282                 do mm=1,ndouble
 283                  ffthet(mm,m,l,i,j,k)=0.0d0
 284                  ggthet(mm,m,l,i,j,k)=0.0d0
 285                 enddo
 286               enddo
 287             enddo
 288           enddo
 289         enddo
 290       enddo
 291       do i=1,nthetyp
 292         do j=1,nthetyp
 293           do k=1,nthetyp
 294             read (ithep,'(3a)') res1,res2,res3
 295             read (ithep,*) aa0thet(i,j,k)
 296             read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
 297             read (ithep,*)
 298      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
 299      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
 300      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
 301      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
 302             read (ithep,*)
 303      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
 304      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
 305      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 306           enddo
 307         enddo
 308       enddo
 309 C
 310 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 311 C coefficients of theta-and-gamma-dependent terms are zero.
 312 C
 313       do i=1,nthetyp
 314         do j=1,nthetyp
 315           do l=1,ntheterm
 316             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
 317             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
 318           enddo
 319           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
 320           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
 321         enddo
 322         do l=1,ntheterm
 323           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
 324         enddo
 325         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
 326       enddo
 327 C
 328 C Control printout of the coefficients of virtual-bond-angle potentials
 329 C
 330       if (lprint) then
 331         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 332         do i=1,nthetyp+1
 333           do j=1,nthetyp+1
 334             do k=1,nthetyp+1
 335               write (iout,'(//4a)')
 336      &         'Type ',onelett(i),onelett(j),onelett(k)
 337               write (iout,'(//a,10x,a)') " l","a[l]"
 338               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
 339               write (iout,'(i2,1pe15.5)')
 340      &           (l,aathet(l,i,j,k),l=1,ntheterm)
 341             do l=1,ntheterm2
 342               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 343      &          "b",l,"c",l,"d",l,"e",l
 344               do m=1,nsingle
 345                 write (iout,'(i2,4(1pe15.5))') m,
 346      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
 347      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
 348               enddo
 349             enddo
 350             do l=1,ntheterm3
 351               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 352      &          "f+",l,"f-",l,"g+",l,"g-",l
 353               do m=2,ndouble
 354                 do n=1,m-1
 355                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 356      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
 357      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
 358                 enddo
 359               enddo
 360             enddo
 361           enddo
 362         enddo
 363       enddo
 364       call flush(iout)
 365       endif
 366 #endif
 367
 368 #ifdef CRYST_SC
 369 C
 370 C Read the parameters of the probability distribution/energy expression
 371 C of the side chains.
 372 C
 373       do i=1,ntyp
 374         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 375         if (i.eq.10) then
 376           dsc_inv(i)=0.0D0
 377         else
 378           dsc_inv(i)=1.0D0/dsc(i)
 379         endif
 380         if (i.ne.10) then
 381         do j=1,nlob(i)
 382           do k=1,3
 383             do l=1,3
 384               blower(l,k,j)=0.0D0
 385             enddo
 386           enddo
 387         enddo
 388         bsc(1,i)=0.0D0
 389         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 390         do j=2,nlob(i)
 391           read (irotam,*) bsc(j,i)
 392           read (irotam,*) (censc(k,j,i),k=1,3),
 393      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 394         enddo
 395         do j=1,nlob(i)
 396           do k=1,3
 397             do l=1,k
 398               akl=0.0D0
 399               do m=1,3
 400                 akl=akl+blower(k,m,j)*blower(l,m,j)
 401               enddo
 402               gaussc(k,l,j,i)=akl
 403               gaussc(l,k,j,i)=akl
 404             enddo
 405           enddo
 406         enddo
 407         endif
 408       enddo
 409       close (irotam)
 410       if (lprint) then
 411         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 412         do i=1,ntyp
 413           nlobi=nlob(i)
 414           if (nlobi.gt.0) then
 415           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 416      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 417 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 418 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 419 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 420 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 421            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 422      &                             'log h',(bsc(j,i),j=1,nlobi)
 423            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 424      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 425 c          write (iout,'(a)')
 426 c         do j=1,nlobi
 427 c           ind=0
 428 c           do k=1,3
 429 c             do l=1,k
 430 c              ind=ind+1
 431 c              blower(k,l,j)=gaussc(ind,j,i)
 432 c             enddo
 433 c           enddo
 434 c         enddo
 435           do k=1,3
 436             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 437      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 438           enddo
 439           endif
 440         enddo
 441       endif
 442 #else
 443 C
 444 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 445 C added by Urszula Kozlowska 07/11/2007
 446 C
 447       do i=1,ntyp
 448         read (irotam,*)
 449        if (i.eq.10) then
 450          read (irotam,*)
 451        else
 452          do j=1,65
 453            read(irotam,*) sc_parmin(j,i)
 454          enddo
 455        endif
 456       enddo
 457 #endif
 458       close(irotam)
 459 #ifdef CRYST_TOR
 460 C
 461 C Read torsional parameters in old format
 462 C
 463       read (itorp,*) ntortyp,nterm_old
 464       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 465       read (itorp,*) (itortyp(i),i=1,ntyp)
 466       do i=1,ntortyp
 467         do j=1,ntortyp
 468           read (itorp,'(a)')
 469           do k=1,nterm_old
 470             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 471           enddo
 472         enddo
 473       enddo
 474       close (itorp)
 475       if (lprint) then
 476         write (iout,'(/a/)') 'Torsional constants:'
 477         do i=1,ntortyp
 478           do j=1,ntortyp
 479             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 480             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 481           enddo
 482         enddo
 483       endif
 484
 485
 486 #else
 487 C
 488 C Read torsional parameters
 489 C
 490       read (itorp,*) ntortyp
 491       read (itorp,*) (itortyp(i),i=1,ntyp)
 492       write (iout,*) 'ntortyp',ntortyp
 493       do i=1,ntortyp
 494         do j=1,ntortyp
 495           read (itorp,*) nterm(i,j),nlor(i,j)
 496           v0ij=0.0d0
 497           si=-1.0d0
 498           do k=1,nterm(i,j)
 499             read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
 500             v0ij=v0ij+si*v1(k,i,j)
 501             si=-si
 502           enddo
 503           do k=1,nlor(i,j)
 504             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 505             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 506           enddo
 507           v0(i,j)=v0ij
 508         enddo
 509       enddo
 510       close (itorp)
 511       if (lprint) then
 512         write (iout,'(/a/)') 'Torsional constants:'
 513         do i=1,ntortyp
 514           do j=1,ntortyp
 515             write (iout,*) 'ityp',i,' jtyp',j
 516             write (iout,*) 'Fourier constants'
 517             do k=1,nterm(i,j)
 518               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 519             enddo
 520             write (iout,*) 'Lorenz constants'
 521             do k=1,nlor(i,j)
 522               write (iout,'(3(1pe15.5))')
 523      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 524             enddo
 525           enddo
 526         enddo
 527       endif
 528 C
 529 C 6/23/01 Read parameters for double torsionals
 530 C
 531       do i=1,ntortyp
 532         do j=1,ntortyp
 533           do k=1,ntortyp
 534             read (itordp,'(3a1)') t1,t2,t3
 535             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 536      &        .or. t3.ne.onelett(k)) then
 537               write (iout,*) "Error in double torsional parameter file",
 538      &         i,j,k,t1,t2,t3
 539                stop "Error in double torsional parameter file"
 540             endif
 541             read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
 542             read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
 543             read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
 544             read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
 545             read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
 546             read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
 547      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
 548           enddo
 549         enddo
 550       enddo
 551       if (lprint) then
 552       write (iout,*)
 553       write (iout,*) 'Constants for double torsionals'
 554       do i=1,ntortyp
 555         do j=1,ntortyp
 556           do k=1,ntortyp
 557             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 558      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 559             write (iout,*)
 560             write (iout,*) 'Single angles:'
 561             do l=1,ntermd_1(i,j,k)
 562               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 563      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 564      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 565             enddo
 566             write (iout,*)
 567             write (iout,*) 'Pairs of angles:'
 568             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 569             do l=1,ntermd_2(i,j,k)
 570               write (iout,'(i5,20f10.5)')
 571      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 572             enddo
 573             write (iout,*)
 574             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 575             do l=1,ntermd_2(i,j,k)
 576               write (iout,'(i5,20f10.5)')
 577      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 578             enddo
 579             write (iout,*)
 580           enddo
 581         enddo
 582       enddo
 583       endif
 584 #endif
 585 C
 586 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 587 C         correlation energies.
 588 C
 589       read (isccor,*) nterm_sccor
 590       do i=1,20
 591         do j=1,20
 592           read (isccor,'(a)')
 593           do k=1,nterm_sccor
 594             read (isccor,*)
 595      &        kk,v1sccor(k,i,j),v2sccor(k,i,j)
 596           enddo
 597         enddo
 598       enddo
 599       close (isccor)
 600       if (lprint) then
 601         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 602         do i=1,20
 603           do j=1,20
 604             write (iout,*) 'ityp',i,' jtyp',j
 605             do k=1,nterm_sccor
 606               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
 607             enddo
 608           enddo
 609         enddo
 610       endif
 611 C
 612 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 613 C         interaction energy of the Gly, Ala, and Pro prototypes.
 614 C
 615       read (ifourier,*) nloctyp
 616       do i=1,nloctyp
 617         read (ifourier,*)
 618         read (ifourier,*) (b(ii,i),ii=1,13)
 619         if (lprint) then
 620         write (iout,*) 'Type',i
 621         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 622         endif
 623         B1(1,i)  = b(3,i)
 624         B1(2,i)  = b(5,i)
 625         B1tilde(1,i) = b(3,i)
 626         B1tilde(2,i) =-b(5,i)
 627         B2(1,i)  = b(2,i)
 628         B2(2,i)  = b(4,i)
 629         CC(1,1,i)= b(7,i)
 630         CC(2,2,i)=-b(7,i)
 631         CC(2,1,i)= b(9,i)
 632         CC(1,2,i)= b(9,i)
 633         Ctilde(1,1,i)=b(7,i)
 634         Ctilde(1,2,i)=b(9,i)
 635         Ctilde(2,1,i)=-b(9,i)
 636         Ctilde(2,2,i)=b(7,i)
 637         DD(1,1,i)= b(6,i)
 638         DD(2,2,i)=-b(6,i)
 639         DD(2,1,i)= b(8,i)
 640         DD(1,2,i)= b(8,i)
 641         Dtilde(1,1,i)=b(6,i)
 642         Dtilde(1,2,i)=b(8,i)
 643         Dtilde(2,1,i)=-b(8,i)
 644         Dtilde(2,2,i)=b(6,i)
 645         EE(1,1,i)= b(10,i)+b(11,i)
 646         EE(2,2,i)=-b(10,i)+b(11,i)
 647         EE(2,1,i)= b(12,i)-b(13,i)
 648         EE(1,2,i)= b(12,i)+b(13,i)
 649       enddo
 650       if (lprint) then
 651       do i=1,nloctyp
 652         write (iout,*) 'Type',i
 653         write (iout,*) 'B1'
 654 c        write (iout,'(f10.5)') B1(:,i)
 655         write(iout,*) B1(1,i),B1(2,i)
 656         write (iout,*) 'B2'
 657 c        write (iout,'(f10.5)') B2(:,i)
 658         write(iout,*) B2(1,i),B2(2,i)
 659         write (iout,*) 'CC'
 660         do j=1,2
 661           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 662         enddo
 663         write(iout,*) 'DD'
 664         do j=1,2
 665           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 666         enddo
 667         write(iout,*) 'EE'
 668         do j=1,2
 669           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 670         enddo
 671       enddo
 672       endif
 673 C
 674 C Read electrostatic-interaction parameters
 675 C
 676       if (lprint) then
 677         write (iout,'(/a)') 'Electrostatic interaction constants:'
 678         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 679      &            'IT','JT','APP','BPP','AEL6','AEL3'
 680       endif
 681       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 682       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 683       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 684       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 685       close (ielep)
 686       do i=1,2
 687         do j=1,2
 688         rri=rpp(i,j)**6
 689         app (i,j)=epp(i,j)*rri*rri
 690         bpp (i,j)=-2.0D0*epp(i,j)*rri
 691         ael6(i,j)=elpp6(i,j)*4.2D0**6
 692         ael3(i,j)=elpp3(i,j)*4.2D0**3
 693         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 694      &                    ael6(i,j),ael3(i,j)
 695         enddo
 696       enddo
 697 C
 698 C Read side-chain interaction parameters.
 699 C
 700       read (isidep,*) ipot,expon
 701       if (ipot.lt.1 .or. ipot.gt.5) then
 702         write (iout,'(2a)') 'Error while reading SC interaction',
 703      &               'potential file - unknown potential type.'
 704         stop
 705       endif
 706       expon2=expon/2
 707       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 708      & ', exponents are ',expon,2*expon
 709       goto (10,20,30,30,40) ipot
 710 C----------------------- LJ potential ---------------------------------
 711    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 712       if (lprint) then
 713         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 714         write (iout,'(a/)') 'The epsilon array:'
 715         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 716         write (iout,'(/a)') 'One-body parameters:'
 717         write (iout,'(a,4x,a)') 'residue','sigma'
 718         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 719       endif
 720       goto 50
 721 C----------------------- LJK potential --------------------------------
 722    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 723      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 724       if (lprint) then
 725         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 726         write (iout,'(a/)') 'The epsilon array:'
 727         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 728         write (iout,'(/a)') 'One-body parameters:'
 729         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 730         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 731      &        i=1,ntyp)
 732       endif
 733       goto 50
 734 C---------------------- GB or BP potential -----------------------------
 735    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 736      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 737      &  (alp(i),i=1,ntyp)
 738 C For the GB potential convert sigma'**2 into chi'
 739       if (ipot.eq.4) then
 740         do i=1,ntyp
 741           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 742         enddo
 743       endif
 744       if (lprint) then
 745         write (iout,'(/a/)') 'Parameters of the BP potential:'
 746         write (iout,'(a/)') 'The epsilon array:'
 747         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 748         write (iout,'(/a)') 'One-body parameters:'
 749         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 750      &       '    chip  ','    alph  '
 751         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 752      &                     chip(i),alp(i),i=1,ntyp)
 753       endif
 754       goto 50
 755 C--------------------- GBV potential -----------------------------------
 756    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 757      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 758      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 759       if (lprint) then
 760         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 761         write (iout,'(a/)') 'The epsilon array:'
 762         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 763         write (iout,'(/a)') 'One-body parameters:'
 764         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 765      &      's||/s_|_^2','    chip  ','    alph  '
 766         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 767      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 768       endif
 769    50 continue
 770       close (isidep)
 771 C-----------------------------------------------------------------------
 772 C Calculate the "working" parameters of SC interactions.
 773       do i=2,ntyp
 774         do j=1,i-1
 775           eps(i,j)=eps(j,i)
 776         enddo
 777       enddo
 778       do i=1,ntyp
 779         do j=i,ntyp
 780           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 781           sigma(j,i)=sigma(i,j)
 782           rs0(i,j)=dwa16*sigma(i,j)
 783           rs0(j,i)=rs0(i,j)
 784         enddo
 785       enddo
 786       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 787      & 'Working parameters of the SC interactions:',
 788      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 789      & '  chi1   ','   chi2   '
 790       do i=1,ntyp
 791         do j=i,ntyp
 792           epsij=eps(i,j)
 793           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 794             rrij=sigma(i,j)
 795           else
 796             rrij=rr0(i)+rr0(j)
 797           endif
 798           r0(i,j)=rrij
 799           r0(j,i)=rrij
 800           rrij=rrij**expon
 801           epsij=eps(i,j)
 802           sigeps=dsign(1.0D0,epsij)
 803           epsij=dabs(epsij)
 804           aa(i,j)=epsij*rrij*rrij
 805           bb(i,j)=-sigeps*epsij*rrij
 806           aa(j,i)=aa(i,j)
 807           bb(j,i)=bb(i,j)
 808           if (ipot.gt.2) then
 809             sigt1sq=sigma0(i)**2
 810             sigt2sq=sigma0(j)**2
 811             sigii1=sigii(i)
 812             sigii2=sigii(j)
 813             ratsig1=sigt2sq/sigt1sq
 814             ratsig2=1.0D0/ratsig1
 815             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 816             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 817             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 818           else
 819             rsum_max=sigma(i,j)
 820           endif
 821 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 822             sigmaii(i,j)=rsum_max
 823             sigmaii(j,i)=rsum_max
 824 c         else
 825 c           sigmaii(i,j)=r0(i,j)
 826 c           sigmaii(j,i)=r0(i,j)
 827 c         endif
 828 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 829           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 830             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 831             augm(i,j)=epsij*r_augm**(2*expon)
 832 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 833             augm(j,i)=augm(i,j)
 834           else
 835             augm(i,j)=0.0D0
 836             augm(j,i)=0.0D0
 837           endif
 838           if (lprint) then
 839             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 840      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 841      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 842           endif
 843         enddo
 844       enddo
 845 C
 846 C Define the SC-p interaction constants
 847 C
 848 #ifdef OLDSCP
 849       do i=1,20
 850 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 851 C helix formation)
 852 c       aad(i,1)=0.3D0*4.0D0**12
 853 C Following line for constants currently implemented
 854 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 855         aad(i,1)=1.5D0*4.0D0**12
 856 c       aad(i,1)=0.17D0*5.6D0**12
 857         aad(i,2)=aad(i,1)
 858 C "Soft" SC-p repulsion
 859         bad(i,1)=0.0D0
 860 C Following line for constants currently implemented
 861 c       aad(i,1)=0.3D0*4.0D0**6
 862 C "Hard" SC-p repulsion
 863         bad(i,1)=3.0D0*4.0D0**6
 864 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 865         bad(i,2)=bad(i,1)
 866 c       aad(i,1)=0.0D0
 867 c       aad(i,2)=0.0D0
 868 c       bad(i,1)=1228.8D0
 869 c       bad(i,2)=1228.8D0
 870       enddo
 871 #else
 872 C
 873 C 8/9/01 Read the SC-p interaction constants from file
 874 C
 875       do i=1,ntyp
 876         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 877       enddo
 878       do i=1,ntyp
 879         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 880         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 881         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 882         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 883       enddo
 884
 885       if (lprint) then
 886         write (iout,*) "Parameters of SC-p interactions:"
 887         do i=1,20
 888           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 889      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 890         enddo
 891       endif
 892 #endif
 893 C
 894 C Define the constants of the disulfide bridge
 895 C
 896       ebr=-5.50D0
 897 c
 898 c Old arbitrary potential - commented out.
 899 c
 900 c      dbr= 4.20D0
 901 c      fbr= 3.30D0
 902 c
 903 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 904 c energy surface of diethyl disulfide.
 905 c A. Liwo and U. Kozlowska, 11/24/03
 906 c
 907       D0CM = 3.78d0
 908       AKCM = 15.1d0
 909       AKTH = 11.0d0
 910       AKCT = 12.0d0
 911       V1SS =-1.08d0
 912       V2SS = 7.61d0
 913       V3SS = 13.7d0
 914
 915       if (lprint) then
 916       write (iout,'(/a)') "Disulfide bridge parameters:"
 917       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 918       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 919       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 920       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 921      & ' v3ss:',v3ss
 922       endif
 923       return
 924       end