source/wham/src/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 Read of Side-chain backbone correlation parameters
 586 C Modified 11 May 2012 by Adasko
 587 CCC
 588 C
 589       read (isccor,*) nsccortyp
 590       read (isccor,*) (isccortyp(i),i=1,ntyp)
 591 c      write (iout,*) 'ntortyp',ntortyp
 592       maxinter=3
 593 cc maxinter is maximum interaction sites
 594       do l=1,maxinter
 595       do i=1,nsccortyp
 596         do j=1,nsccortyp
 597           read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
 598           v0ijsccor=0.0d0
 599           si=-1.0d0
 600
 601           do k=1,nterm_sccor(i,j)
 602             read (isccor,*) kk,v1sccor(k,l,i,j)
 603      &    ,v2sccor(k,l,i,j)
 604             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 605             si=-si
 606           enddo
 607           do k=1,nlor_sccor(i,j)
 608             read (isccor,*) kk,vlor1sccor(k,i,j),
 609      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 610             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 611      &(1+vlor3sccor(k,i,j)**2)
 612           enddo
 613           v0sccor(i,j)=v0ijsccor
 614         enddo
 615       enddo
 616       enddo
 617       close (isccor)
 618
 619       if (lprint) then
 620         write (iout,'(/a/)') 'Torsional constants:'
 621         do i=1,nsccortyp
 622           do j=1,nsccortyp
 623             write (iout,*) 'ityp',i,' jtyp',j
 624             write (iout,*) 'Fourier constants'
 625             do k=1,nterm_sccor(i,j)
 626            write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 627             enddo
 628             write (iout,*) 'Lorenz constants'
 629             do k=1,nlor_sccor(i,j)
 630               write (iout,'(3(1pe15.5))')
 631      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 632             enddo
 633           enddo
 634         enddo
 635       endif
 636
 637 C
 638 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 639 C         interaction energy of the Gly, Ala, and Pro prototypes.
 640 C
 641       read (ifourier,*) nloctyp
 642       do i=1,nloctyp
 643         read (ifourier,*)
 644         read (ifourier,*) (b(ii,i),ii=1,13)
 645         if (lprint) then
 646         write (iout,*) 'Type',i
 647         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 648         endif
 649         B1(1,i)  = b(3,i)
 650         B1(2,i)  = b(5,i)
 651         B1tilde(1,i) = b(3,i)
 652         B1tilde(2,i) =-b(5,i)
 653         B2(1,i)  = b(2,i)
 654         B2(2,i)  = b(4,i)
 655         CC(1,1,i)= b(7,i)
 656         CC(2,2,i)=-b(7,i)
 657         CC(2,1,i)= b(9,i)
 658         CC(1,2,i)= b(9,i)
 659         Ctilde(1,1,i)=b(7,i)
 660         Ctilde(1,2,i)=b(9,i)
 661         Ctilde(2,1,i)=-b(9,i)
 662         Ctilde(2,2,i)=b(7,i)
 663         DD(1,1,i)= b(6,i)
 664         DD(2,2,i)=-b(6,i)
 665         DD(2,1,i)= b(8,i)
 666         DD(1,2,i)= b(8,i)
 667         Dtilde(1,1,i)=b(6,i)
 668         Dtilde(1,2,i)=b(8,i)
 669         Dtilde(2,1,i)=-b(8,i)
 670         Dtilde(2,2,i)=b(6,i)
 671         EE(1,1,i)= b(10,i)+b(11,i)
 672         EE(2,2,i)=-b(10,i)+b(11,i)
 673         EE(2,1,i)= b(12,i)-b(13,i)
 674         EE(1,2,i)= b(12,i)+b(13,i)
 675       enddo
 676       if (lprint) then
 677       do i=1,nloctyp
 678         write (iout,*) 'Type',i
 679         write (iout,*) 'B1'
 680 c        write (iout,'(f10.5)') B1(:,i)
 681         write(iout,*) B1(1,i),B1(2,i)
 682         write (iout,*) 'B2'
 683 c        write (iout,'(f10.5)') B2(:,i)
 684         write(iout,*) B2(1,i),B2(2,i)
 685         write (iout,*) 'CC'
 686         do j=1,2
 687           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 688         enddo
 689         write(iout,*) 'DD'
 690         do j=1,2
 691           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 692         enddo
 693         write(iout,*) 'EE'
 694         do j=1,2
 695           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 696         enddo
 697       enddo
 698       endif
 699 C
 700 C Read electrostatic-interaction parameters
 701 C
 702       if (lprint) then
 703         write (iout,'(/a)') 'Electrostatic interaction constants:'
 704         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 705      &            'IT','JT','APP','BPP','AEL6','AEL3'
 706       endif
 707       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 708       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 709       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 710       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 711       close (ielep)
 712       do i=1,2
 713         do j=1,2
 714         rri=rpp(i,j)**6
 715         app (i,j)=epp(i,j)*rri*rri
 716         bpp (i,j)=-2.0D0*epp(i,j)*rri
 717         ael6(i,j)=elpp6(i,j)*4.2D0**6
 718         ael3(i,j)=elpp3(i,j)*4.2D0**3
 719         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 720      &                    ael6(i,j),ael3(i,j)
 721         enddo
 722       enddo
 723 C
 724 C Read side-chain interaction parameters.
 725 C
 726       read (isidep,*) ipot,expon
 727       if (ipot.lt.1 .or. ipot.gt.5) then
 728         write (iout,'(2a)') 'Error while reading SC interaction',
 729      &               'potential file - unknown potential type.'
 730         stop
 731       endif
 732       expon2=expon/2
 733       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 734      & ', exponents are ',expon,2*expon
 735       goto (10,20,30,30,40) ipot
 736 C----------------------- LJ potential ---------------------------------
 737    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 738       if (lprint) then
 739         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 740         write (iout,'(a/)') 'The epsilon array:'
 741         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 742         write (iout,'(/a)') 'One-body parameters:'
 743         write (iout,'(a,4x,a)') 'residue','sigma'
 744         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 745       endif
 746       goto 50
 747 C----------------------- LJK potential --------------------------------
 748    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 749      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 750       if (lprint) then
 751         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 752         write (iout,'(a/)') 'The epsilon array:'
 753         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 754         write (iout,'(/a)') 'One-body parameters:'
 755         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 756         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 757      &        i=1,ntyp)
 758       endif
 759       goto 50
 760 C---------------------- GB or BP potential -----------------------------
 761    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 762      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 763      &  (alp(i),i=1,ntyp)
 764 C For the GB potential convert sigma'**2 into chi'
 765       if (ipot.eq.4) then
 766         do i=1,ntyp
 767           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 768         enddo
 769       endif
 770       if (lprint) then
 771         write (iout,'(/a/)') 'Parameters of the BP potential:'
 772         write (iout,'(a/)') 'The epsilon array:'
 773         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 774         write (iout,'(/a)') 'One-body parameters:'
 775         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 776      &       '    chip  ','    alph  '
 777         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 778      &                     chip(i),alp(i),i=1,ntyp)
 779       endif
 780       goto 50
 781 C--------------------- GBV potential -----------------------------------
 782    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 783      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 784      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 785       if (lprint) then
 786         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 787         write (iout,'(a/)') 'The epsilon array:'
 788         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 789         write (iout,'(/a)') 'One-body parameters:'
 790         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 791      &      's||/s_|_^2','    chip  ','    alph  '
 792         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 793      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 794       endif
 795    50 continue
 796       close (isidep)
 797 C-----------------------------------------------------------------------
 798 C Calculate the "working" parameters of SC interactions.
 799       do i=2,ntyp
 800         do j=1,i-1
 801           eps(i,j)=eps(j,i)
 802         enddo
 803       enddo
 804       do i=1,ntyp
 805         do j=i,ntyp
 806           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 807           sigma(j,i)=sigma(i,j)
 808           rs0(i,j)=dwa16*sigma(i,j)
 809           rs0(j,i)=rs0(i,j)
 810         enddo
 811       enddo
 812       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 813      & 'Working parameters of the SC interactions:',
 814      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 815      & '  chi1   ','   chi2   '
 816       do i=1,ntyp
 817         do j=i,ntyp
 818           epsij=eps(i,j)
 819           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 820             rrij=sigma(i,j)
 821           else
 822             rrij=rr0(i)+rr0(j)
 823           endif
 824           r0(i,j)=rrij
 825           r0(j,i)=rrij
 826           rrij=rrij**expon
 827           epsij=eps(i,j)
 828           sigeps=dsign(1.0D0,epsij)
 829           epsij=dabs(epsij)
 830           aa(i,j)=epsij*rrij*rrij
 831           bb(i,j)=-sigeps*epsij*rrij
 832           aa(j,i)=aa(i,j)
 833           bb(j,i)=bb(i,j)
 834           if (ipot.gt.2) then
 835             sigt1sq=sigma0(i)**2
 836             sigt2sq=sigma0(j)**2
 837             sigii1=sigii(i)
 838             sigii2=sigii(j)
 839             ratsig1=sigt2sq/sigt1sq
 840             ratsig2=1.0D0/ratsig1
 841             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 842             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 843             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 844           else
 845             rsum_max=sigma(i,j)
 846           endif
 847 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 848             sigmaii(i,j)=rsum_max
 849             sigmaii(j,i)=rsum_max
 850 c         else
 851 c           sigmaii(i,j)=r0(i,j)
 852 c           sigmaii(j,i)=r0(i,j)
 853 c         endif
 854 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 855           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 856             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 857             augm(i,j)=epsij*r_augm**(2*expon)
 858 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 859             augm(j,i)=augm(i,j)
 860           else
 861             augm(i,j)=0.0D0
 862             augm(j,i)=0.0D0
 863           endif
 864           if (lprint) then
 865             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 866      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 867      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 868           endif
 869         enddo
 870       enddo
 871 C
 872 C Define the SC-p interaction constants
 873 C
 874 #ifdef OLDSCP
 875       do i=1,20
 876 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 877 C helix formation)
 878 c       aad(i,1)=0.3D0*4.0D0**12
 879 C Following line for constants currently implemented
 880 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 881         aad(i,1)=1.5D0*4.0D0**12
 882 c       aad(i,1)=0.17D0*5.6D0**12
 883         aad(i,2)=aad(i,1)
 884 C "Soft" SC-p repulsion
 885         bad(i,1)=0.0D0
 886 C Following line for constants currently implemented
 887 c       aad(i,1)=0.3D0*4.0D0**6
 888 C "Hard" SC-p repulsion
 889         bad(i,1)=3.0D0*4.0D0**6
 890 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 891         bad(i,2)=bad(i,1)
 892 c       aad(i,1)=0.0D0
 893 c       aad(i,2)=0.0D0
 894 c       bad(i,1)=1228.8D0
 895 c       bad(i,2)=1228.8D0
 896       enddo
 897 #else
 898 C
 899 C 8/9/01 Read the SC-p interaction constants from file
 900 C
 901       do i=1,ntyp
 902         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 903       enddo
 904       do i=1,ntyp
 905         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 906         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 907         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 908         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 909       enddo
 910
 911       if (lprint) then
 912         write (iout,*) "Parameters of SC-p interactions:"
 913         do i=1,20
 914           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 915      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 916         enddo
 917       endif
 918 #endif
 919 C
 920 C Define the constants of the disulfide bridge
 921 C
 922       ebr=-5.50D0
 923 c
 924 c Old arbitrary potential - commented out.
 925 c
 926 c      dbr= 4.20D0
 927 c      fbr= 3.30D0
 928 c
 929 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 930 c energy surface of diethyl disulfide.
 931 c A. Liwo and U. Kozlowska, 11/24/03
 932 c
 933       D0CM = 3.78d0
 934       AKCM = 15.1d0
 935       AKTH = 11.0d0
 936       AKCT = 12.0d0
 937       V1SS =-1.08d0
 938       V2SS = 7.61d0
 939       V3SS = 13.7d0
 940
 941       if (lprint) then
 942       write (iout,'(/a)') "Disulfide bridge parameters:"
 943       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 944       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 945       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 946       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 947      & ' v3ss:',v3ss
 948       endif
 949       return
 950       end