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