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