source/wham/src-NEWSC-NEWCORR/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=.true.
  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,*,end=110,err=110) vbldp0,akp
 154       do i=1,ntyp
 155         nbondterm(i)=1
 156         read (ibond,*,end=110,err=110) 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,*,end=110,err=110) ijunk,vbldp0,akp,rjunk
 166       do i=1,ntyp
 167         read (ibond,*,end=110,err=110) nbondterm(i),(vbldsc0(j,i),
 168      &   aksc(j,i),abond0(j,i),
 169      &   j=1,nbondterm(i))
 170         dsc(i) = vbldsc0(1,i)
 171         if (i.eq.10) then
 172           dsc_inv(i)=0.0D0
 173         else
 174           dsc_inv(i)=1.0D0/dsc(i)
 175         endif
 176       enddo
 177 #endif
 178       if (lprint) then
 179         write(iout,'(/a/)')"Force constants virtual bonds:"
 180         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
 181      &   'inertia','Pstok'
 182         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 183         do i=1,ntyp
 184           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
 185      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
 186           do j=2,nbondterm(i)
 187             write (iout,'(13x,3f10.5)')
 188      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
 189           enddo
 190         enddo
 191       endif
 192 #ifdef CRYST_THETA
 193 C
 194 C Read the parameters of the probability distribution/energy expression
 195 C of the virtual-bond valence angles theta
 196 C
 197       do i=1,ntyp
 198         read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i),j=1,2),
 199      &     (bthet(j,i),j=1,2)
 200         read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
 201         read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
 202         read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
 203         sigc0(i)=sigc0(i)**2
 204       enddo
 205       close (ithep)
 206       if (lprint) then
 207 c       write (iout,'(a)')
 208 c    &    'Parameters of the virtual-bond valence angles:'
 209 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 210 c    & '    ATHETA0   ','         A1   ','        A2    ',
 211 c    & '        B1    ','         B2   '
 212 c       do i=1,ntyp
 213 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 214 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 215 c       enddo
 216 c       write (iout,'(/a/9x,5a/79(1h-))')
 217 c    & 'Parameters of the expression for sigma(theta_c):',
 218 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 219 c    & '     ALPH3    ','    SIGMA0C   '
 220 c       do i=1,ntyp
 221 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 222 c    &      (polthet(j,i),j=0,3),sigc0(i)
 223 c       enddo
 224 c       write (iout,'(/a/9x,5a/79(1h-))')
 225 c    & 'Parameters of the second gaussian:',
 226 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 227 c    & '        G2    ','         G3   '
 228 c       do i=1,ntyp
 229 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 230 c    &       sig0(i),(gthet(j,i),j=1,3)
 231 c       enddo
 232         write (iout,'(a)')
 233      &    'Parameters of the virtual-bond valence angles:'
 234         write (iout,'(/a/9x,5a/79(1h-))')
 235      & 'Coefficients of expansion',
 236      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 237      & '   b1*10^1    ','    b2*10^1   '
 238         do i=1,ntyp
 239           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 240      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
 241         enddo
 242         write (iout,'(/a/9x,5a/79(1h-))')
 243      & 'Parameters of the expression for sigma(theta_c):',
 244      & ' alpha0       ','  alph1       ',' alph2        ',
 245      & ' alhp3        ','   sigma0c    '
 246         do i=1,ntyp
 247           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 248      &      (polthet(j,i),j=0,3),sigc0(i)
 249         enddo
 250         write (iout,'(/a/9x,5a/79(1h-))')
 251      & 'Parameters of the second gaussian:',
 252      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 253      & '        G2    ','   G3*10^1    '
 254         do i=1,ntyp
 255           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 256      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 257         enddo
 258       endif
 259 #else
 260 C
 261 C Read the parameters of Utheta determined from ab initio surfaces
 262 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 263 C
 264       read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
 265      &  ntheterm3,nsingle,ndouble
 266       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 267       read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
 268       do i=1,maxthetyp
 269         do j=1,maxthetyp
 270           do k=1,maxthetyp
 271             aa0thet(i,j,k)=0.0d0
 272             do l=1,ntheterm
 273               aathet(l,i,j,k)=0.0d0
 274             enddo
 275             do l=1,ntheterm2
 276               do m=1,nsingle
 277                 bbthet(m,l,i,j,k)=0.0d0
 278                 ccthet(m,l,i,j,k)=0.0d0
 279                 ddthet(m,l,i,j,k)=0.0d0
 280                 eethet(m,l,i,j,k)=0.0d0
 281               enddo
 282             enddo
 283             do l=1,ntheterm3
 284               do m=1,ndouble
 285                 do mm=1,ndouble
 286                  ffthet(mm,m,l,i,j,k)=0.0d0
 287                  ggthet(mm,m,l,i,j,k)=0.0d0
 288                 enddo
 289               enddo
 290             enddo
 291           enddo
 292         enddo
 293       enddo
 294       do i=1,nthetyp
 295         do j=1,nthetyp
 296           do k=1,nthetyp
 297             read (ithep,'(3a)',end=111,err=111) res1,res2,res3
 298             read (ithep,*,end=111,err=111) aa0thet(i,j,k)
 299             read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
 300             read (ithep,*,end=111,err=111)
 301      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
 302      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
 303      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
 304      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
 305             read (ithep,*,end=111,err=111)
 306      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
 307      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
 308      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 309           enddo
 310         enddo
 311       enddo
 312 C
 313 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 314 C coefficients of theta-and-gamma-dependent terms are zero.
 315 C
 316       do i=1,nthetyp
 317         do j=1,nthetyp
 318           do l=1,ntheterm
 319             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
 320             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
 321           enddo
 322           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
 323           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
 324         enddo
 325         do l=1,ntheterm
 326           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
 327         enddo
 328         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
 329       enddo
 330 C
 331 C Control printout of the coefficients of virtual-bond-angle potentials
 332 C
 333       if (lprint) then
 334         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 335         do i=1,nthetyp+1
 336           do j=1,nthetyp+1
 337             do k=1,nthetyp+1
 338               write (iout,'(//4a)')
 339      &         'Type ',onelett(i),onelett(j),onelett(k)
 340               write (iout,'(//a,10x,a)') " l","a[l]"
 341               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
 342               write (iout,'(i2,1pe15.5)')
 343      &           (l,aathet(l,i,j,k),l=1,ntheterm)
 344             do l=1,ntheterm2
 345               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 346      &          "b",l,"c",l,"d",l,"e",l
 347               do m=1,nsingle
 348                 write (iout,'(i2,4(1pe15.5))') m,
 349      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
 350      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
 351               enddo
 352             enddo
 353             do l=1,ntheterm3
 354               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 355      &          "f+",l,"f-",l,"g+",l,"g-",l
 356               do m=2,ndouble
 357                 do n=1,m-1
 358                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 359      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
 360      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
 361                 enddo
 362               enddo
 363             enddo
 364           enddo
 365         enddo
 366       enddo
 367       call flush(iout)
 368       endif
 369 #endif
 370
 371 #ifdef CRYST_SC
 372 C
 373 C Read the parameters of the probability distribution/energy expression
 374 C of the side chains.
 375 C
 376       do i=1,ntyp
 377         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 378         if (i.eq.10) then
 379           dsc_inv(i)=0.0D0
 380         else
 381           dsc_inv(i)=1.0D0/dsc(i)
 382         endif
 383         if (i.ne.10) then
 384         do j=1,nlob(i)
 385           do k=1,3
 386             do l=1,3
 387               blower(l,k,j)=0.0D0
 388             enddo
 389           enddo
 390         enddo
 391         bsc(1,i)=0.0D0
 392         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 393      &     ((blower(k,l,1),l=1,k),k=1,3)
 394         do j=2,nlob(i)
 395           read (irotam,*,end=112,err=112) bsc(j,i)
 396           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 397      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 398         enddo
 399         do j=1,nlob(i)
 400           do k=1,3
 401             do l=1,k
 402               akl=0.0D0
 403               do m=1,3
 404                 akl=akl+blower(k,m,j)*blower(l,m,j)
 405               enddo
 406               gaussc(k,l,j,i)=akl
 407               gaussc(l,k,j,i)=akl
 408             enddo
 409           enddo
 410         enddo
 411         endif
 412       enddo
 413       close (irotam)
 414       if (lprint) then
 415         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 416         do i=1,ntyp
 417           nlobi=nlob(i)
 418           if (nlobi.gt.0) then
 419           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 420      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 421 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 422 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 423 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 424 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 425            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 426      &                             'log h',(bsc(j,i),j=1,nlobi)
 427            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 428      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 429 c          write (iout,'(a)')
 430 c         do j=1,nlobi
 431 c           ind=0
 432 c           do k=1,3
 433 c             do l=1,k
 434 c              ind=ind+1
 435 c              blower(k,l,j)=gaussc(ind,j,i)
 436 c             enddo
 437 c           enddo
 438 c         enddo
 439           do k=1,3
 440             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 441      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 442           enddo
 443           endif
 444         enddo
 445       endif
 446 #else
 447 C
 448 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 449 C added by Urszula Kozlowska 07/11/2007
 450 C
 451       do i=1,ntyp
 452         read (irotam,*,end=112,err=112)
 453        if (i.eq.10) then
 454          read (irotam,*,end=112,err=112)
 455        else
 456          do j=1,65
 457            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 458          enddo
 459        endif
 460       enddo
 461 #endif
 462       close(irotam)
 463 #ifdef CRYST_TOR
 464 C
 465 C Read torsional parameters in old format
 466 C
 467       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 468       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 469       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 470       do i=1,ntortyp
 471         do j=1,ntortyp
 472           read (itorp,'(a)',end=113,err=113)
 473           do k=1,nterm_old
 474             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 475           enddo
 476         enddo
 477       enddo
 478       close (itorp)
 479       if (lprint) then
 480         write (iout,'(/a/)') 'Torsional constants:'
 481         do i=1,ntortyp
 482           do j=1,ntortyp
 483             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 484             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 485           enddo
 486         enddo
 487       endif
 488
 489
 490 #else
 491 C
 492 C Read torsional parameters
 493 C
 494       read (itorp,*,end=113,err=113) ntortyp
 495       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 496       write (iout,*) 'ntortyp',ntortyp
 497       do i=1,ntortyp
 498         do j=1,ntortyp
 499           read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
 500           v0ij=0.0d0
 501           si=-1.0d0
 502           do k=1,nterm(i,j)
 503             read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
 504             v0ij=v0ij+si*v1(k,i,j)
 505             si=-si
 506           enddo
 507           do k=1,nlor(i,j)
 508             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),vlor2(k,i,j),
 509      &         vlor3(k,i,j)
 510             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 511           enddo
 512           v0(i,j)=v0ij
 513         enddo
 514       enddo
 515       close (itorp)
 516       if (lprint) then
 517         write (iout,'(/a/)') 'Torsional constants:'
 518         do i=1,ntortyp
 519           do j=1,ntortyp
 520             write (iout,*) 'ityp',i,' jtyp',j
 521             write (iout,*) 'Fourier constants'
 522             do k=1,nterm(i,j)
 523               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 524             enddo
 525             write (iout,*) 'Lorenz constants'
 526             do k=1,nlor(i,j)
 527               write (iout,'(3(1pe15.5))')
 528      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 529             enddo
 530           enddo
 531         enddo
 532       endif
 533 C
 534 C 6/23/01 Read parameters for double torsionals
 535 C
 536       do i=1,ntortyp
 537         do j=1,ntortyp
 538           do k=1,ntortyp
 539             read (itordp,'(3a1)',end=112,err=112) t1,t2,t3
 540             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 541      &        .or. t3.ne.onelett(k)) then
 542               write (iout,*) "Error in double torsional parameter file",
 543      &         i,j,k,t1,t2,t3
 544                stop "Error in double torsional parameter file"
 545             endif
 546             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
 547      &        ntermd_2(i,j,k)
 548             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),
 549      &        l=1,ntermd_1(i,j,k))
 550             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
 551      &        ntermd_1(i,j,k))
 552             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),
 553      &        l=1,ntermd_1(i,j,k))
 554             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
 555      &        ntermd_1(i,j,k))
 556             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
 557      &        v2c(m,l,i,j,k),
 558      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
 559           enddo
 560         enddo
 561       enddo
 562       if (lprint) then
 563       write (iout,*)
 564       write (iout,*) 'Constants for double torsionals'
 565       do i=1,ntortyp
 566         do j=1,ntortyp
 567           do k=1,ntortyp
 568             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 569      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 570             write (iout,*)
 571             write (iout,*) 'Single angles:'
 572             do l=1,ntermd_1(i,j,k)
 573               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 574      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 575      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 576             enddo
 577             write (iout,*)
 578             write (iout,*) 'Pairs of angles:'
 579             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 580             do l=1,ntermd_2(i,j,k)
 581               write (iout,'(i5,20f10.5)')
 582      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 583             enddo
 584             write (iout,*)
 585             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 586             do l=1,ntermd_2(i,j,k)
 587               write (iout,'(i5,20f10.5)')
 588      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 589             enddo
 590             write (iout,*)
 591           enddo
 592         enddo
 593       enddo
 594       endif
 595 #endif
 596 C Read of Side-chain backbone correlation parameters
 597 C Modified 11 May 2012 by Adasko
 598 CCC
 599 C
 600       read (isccor,*,end=119,err=119) nsccortyp
 601       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 602 c      write (iout,*) 'ntortyp',ntortyp
 603       maxinter=3
 604 cc maxinter is maximum interaction sites
 605       do l=1,maxinter
 606       do i=1,nsccortyp
 607         do j=1,nsccortyp
 608           read (isccor,*,end=119,err=119) nterm_sccor(i,j),
 609      &       nlor_sccor(i,j)
 610           v0ijsccor=0.0d0
 611           si=-1.0d0
 612           do k=1,nterm_sccor(i,j)
 613             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 614      &    ,v2sccor(k,l,i,j)
 615             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 616             si=-si
 617           enddo
 618           do k=1,nlor_sccor(i,j)
 619             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 620      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 621             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 622      &(1+vlor3sccor(k,i,j)**2)
 623           enddo
 624           v0sccor(i,j)=v0ijsccor
 625         enddo
 626       enddo
 627       enddo
 628       close (isccor)
 629
 630       if (lprint) then
 631         write (iout,'(/a/)') 'Torsional constants:'
 632         do i=1,nsccortyp
 633           do j=1,nsccortyp
 634             write (iout,*) 'ityp',i,' jtyp',j
 635             write (iout,*) 'Fourier constants'
 636             do k=1,nterm_sccor(i,j)
 637            write(iout,'(2(1pe15.5))')(v1sccor(k,l,i,j),v2sccor(k,l,i,j),
 638      &       l=1,maxinter)
 639             enddo
 640             write (iout,*) 'Lorenz constants'
 641             do k=1,nlor_sccor(i,j)
 642               write (iout,'(3(1pe15.5))')
 643      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 644             enddo
 645           enddo
 646         enddo
 647       endif
 648
 649 C
 650 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 651 C         interaction energy of the Gly, Ala, and Pro prototypes.
 652 C
 653       read (ifourier,*,end=115,err=115) nloctyp
 654       do i=1,nloctyp
 655         read (ifourier,*,end=115,err=115)
 656         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
 657 #ifdef NEWCORR
 658         read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
 659         read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
 660         read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
 661         read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
 662         read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
 663 #endif
 664         if (lprint) then
 665         write (iout,*) 'Type',i
 666         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 667         endif
 668 #ifdef NEWCORR
 669         B1(1,i)  = b(3,i)
 670         B1(2,i)  = b(5,i)
 671         B1tilde(1,i) = b(3,i)
 672         B1tilde(2,i) =-b(5,i)
 673         B2(1,i)  = b(2,i)
 674         B2(2,i)  = b(4,i)
 675 #endif
 676         CC(1,1,i)= b(7,i)
 677         CC(2,2,i)=-b(7,i)
 678         CC(2,1,i)= b(9,i)
 679         CC(1,2,i)= b(9,i)
 680         Ctilde(1,1,i)=b(7,i)
 681         Ctilde(1,2,i)=b(9,i)
 682         Ctilde(2,1,i)=-b(9,i)
 683         Ctilde(2,2,i)=b(7,i)
 684         DD(1,1,i)= b(6,i)
 685         DD(2,2,i)=-b(6,i)
 686         DD(2,1,i)= b(8,i)
 687         DD(1,2,i)= b(8,i)
 688         Dtilde(1,1,i)=b(6,i)
 689         Dtilde(1,2,i)=b(8,i)
 690         Dtilde(2,1,i)=-b(8,i)
 691         Dtilde(2,2,i)=b(6,i)
 692 #ifdef NEWCORR
 693         EEold(1,1,i)= b(10,i)+b(11,i)
 694         EEold(2,2,i)=-b(10,i)+b(11,i)
 695         EEold(2,1,i)= b(12,i)-b(13,i)
 696         EEold(1,2,i)= b(12,i)+b(13,i)
 697         EEold(1,1,-i)= b(10,i)+b(11,i)
 698         EEold(2,2,-i)=-b(10,i)+b(11,i)
 699         EEold(2,1,-i)=-b(12,i)+b(13,i)
 700         EEold(1,2,-i)=-b(12,i)-b(13,i)
 701 #else
 702         EE(1,1,i)= b(10,i)+b(11,i)
 703         EE(2,2,i)=-b(10,i)+b(11,i)
 704         EE(2,1,i)= b(12,i)-b(13,i)
 705         EE(1,2,i)= b(12,i)+b(13,i)
 706 #endif
 707       enddo
 708       if (lprint) then
 709       do i=1,nloctyp
 710         write (iout,*) 'Type',i
 711         write (iout,*) 'B1'
 712 c        write (iout,'(f10.5)') B1(:,i)
 713         write(iout,*) B1(1,i),B1(2,i)
 714         write (iout,*) 'B2'
 715 c        write (iout,'(f10.5)') B2(:,i)
 716         write(iout,*) B2(1,i),B2(2,i)
 717         write (iout,*) 'CC'
 718         do j=1,2
 719           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 720         enddo
 721         write(iout,*) 'DD'
 722         do j=1,2
 723           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 724         enddo
 725         write(iout,*) 'EE'
 726         do j=1,2
 727           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 728         enddo
 729       enddo
 730       endif
 731 C
 732 C Read electrostatic-interaction parameters
 733 C
 734       if (lprint) then
 735         write (iout,'(/a)') 'Electrostatic interaction constants:'
 736         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 737      &            'IT','JT','APP','BPP','AEL6','AEL3'
 738       endif
 739       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 740       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 741       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 742       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 743       close (ielep)
 744       do i=1,2
 745         do j=1,2
 746         rri=rpp(i,j)**6
 747         app (i,j)=epp(i,j)*rri*rri
 748         bpp (i,j)=-2.0D0*epp(i,j)*rri
 749         ael6(i,j)=elpp6(i,j)*4.2D0**6
 750         ael3(i,j)=elpp3(i,j)*4.2D0**3
 751         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 752      &                    ael6(i,j),ael3(i,j)
 753         enddo
 754       enddo
 755 C
 756 C Read side-chain interaction parameters.
 757 C
 758       read (isidep,*,end=117,err=117) ipot,expon
 759       if (ipot.lt.1 .or. ipot.gt.6) then
 760         write (iout,'(2a)') 'Error while reading SC interaction',
 761      &               'potential file - unknown potential type.'
 762         stop
 763       endif
 764       expon2=expon/2
 765       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 766      & ', exponents are ',expon,2*expon
 767       goto (10,20,30,30,40,50) ipot
 768 C----------------------- LJ potential ---------------------------------
 769    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
 770      &        (sigma0(i),i=1,ntyp)
 771       if (lprint) then
 772         write (iout,'(/a/)') 'Parameters of the LJ 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,a)') 'residue','sigma'
 777         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 778       endif
 779       goto 60
 780 C----------------------- LJK potential --------------------------------
 781    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 782      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 783       if (lprint) then
 784         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 785         write (iout,'(a/)') 'The epsilon array:'
 786         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 787         write (iout,'(/a)') 'One-body parameters:'
 788         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 789         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 790      &        i=1,ntyp)
 791       endif
 792       goto 60
 793 C---------------------- GB or BP potential -----------------------------
 794    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 795      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 796      &  (alp(i),i=1,ntyp)
 797 C For the GB potential convert sigma'**2 into chi'
 798       if (ipot.eq.4) then
 799         do i=1,ntyp
 800           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 801         enddo
 802       endif
 803       if (lprint) then
 804         write (iout,'(/a/)') 'Parameters of the BP potential:'
 805         write (iout,'(a/)') 'The epsilon array:'
 806         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 807         write (iout,'(/a)') 'One-body parameters:'
 808         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 809      &       '    chip  ','    alph  '
 810         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 811      &                     chip(i),alp(i),i=1,ntyp)
 812       endif
 813       goto 60
 814 C--------------------- GBV potential -----------------------------------
 815    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 816      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 817      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 818       if (lprint) then
 819         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 820         write (iout,'(a/)') 'The epsilon array:'
 821         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 822         write (iout,'(/a)') 'One-body parameters:'
 823         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 824      &      's||/s_|_^2','    chip  ','    alph  '
 825         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 826      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 827       endif
 828       goto 60
 829 C--------------------- Momo potential -----------------------------------
 830
 831    50 continue
 832
 833       read (isidep,*,end=116,err=116) (icharge(i),i=1,ntyp)
 834 c      write (2,*) "icharge",(icharge(i),i=1,ntyp)
 835       do i=1,ntyp
 836        do j=1,i
 837 c!        write (*,*) "Im in ", i, " ", j
 838         read(isidep,*,end=116,err=116)
 839      &  eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
 840      &  (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
 841      &  chis(i,j),chis(j,i),
 842      &  nstate(i,j),(wstate(k,i,j),k=1,4),
 843      &  dhead(1,1,i,j),
 844      &  dhead(1,2,i,j),
 845      &  dhead(2,1,i,j),
 846      &  dhead(2,2,i,j),
 847      &  dtail(1,i,j),dtail(2,i,j),
 848      &  epshead(i,j),sig0head(i,j),
 849      &  rborn(i,j),rborn(j,i),
 850      &  (wqdip(k,i,j),k=1,2),wquad(i,j),
 851      &  alphapol(i,j),alphapol(j,i),
 852      &  (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
 853        END DO
 854       END DO
 855 c!      write (*,*) "nstate gly-gly", nstate(10,10)
 856 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
 857 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
 858 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
 859 c! FROM HIS-ARG.
 860 c!
 861 c! DISABLE IT AT >>YOUR OWN PERIL<<
 862 c!
 863       DO i = 1, ntyp
 864        DO j = i+1, ntyp
 865         eps(i,j) = eps(j,i)
 866         sigma(i,j) = sigma(j,i)
 867         nstate(i,j) = nstate(j,i)
 868         sigmap1(i,j) = sigmap1(j,i)
 869         sigmap2(i,j) = sigmap2(j,i)
 870         sigiso1(i,j) = sigiso1(j,i)
 871         sigiso2(i,j) = sigiso2(j,i)
 872
 873         DO k = 1, 4
 874          alphasur(k,i,j) = alphasur(k,j,i)
 875          wstate(k,i,j) = wstate(k,j,i)
 876          alphiso(k,i,j) = alphiso(k,j,i)
 877         END DO
 878
 879         dhead(2,1,i,j) = dhead(1,1,j,i)
 880         dhead(2,2,i,j) = dhead(1,2,j,i)
 881         dhead(1,1,i,j) = dhead(2,1,j,i)
 882         dhead(1,2,i,j) = dhead(2,2,j,i)
 883         dtail(1,i,j) = dtail(1,j,i)
 884         dtail(2,i,j) = dtail(2,j,i)
 885 c!        DO k = 1, 2
 886 c!         DO l = 1, 2
 887 c!         write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
 888 c!         write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
 889 c!          dhead(l,k,i,j) = dhead(k,l,j,i)
 890 c!         END DO
 891 c!        END DO
 892
 893         epshead(i,j) = epshead(j,i)
 894         sig0head(i,j) = sig0head(j,i)
 895
 896         DO k = 1, 2
 897          wqdip(k,i,j) = wqdip(k,j,i)
 898         END DO
 899
 900         wquad(i,j) = wquad(j,i)
 901         epsintab(i,j) = epsintab(j,i)
 902
 903        END DO
 904       END DO
 905
 906       if (.not.lprint) goto 70
 907       write (iout,'(a)')
 908      & "Parameters of the new physics-based SC-SC interaction potential"
 909       write (iout,'(/7a)') 'Residues','     epsGB','       rGB',
 910      &  '    chi1GB','    chi2GB','   chip1GB','   chip2GB'
 911       do i=1,ntyp
 912         do j=1,i
 913           write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
 914      &      restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
 915      &       chipp(i,j),chipp(j,i)
 916         enddo
 917       enddo
 918       write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
 919      &  ' alphasur3',' alphasur4','   sigmap1','   sigmap2',
 920      &  '     chis1','     chis2'
 921       do i=1,ntyp
 922         do j=1,i
 923           write (iout,'(2(a3,1x),8(0pf10.3))')
 924      &      restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
 925      &       sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
 926         enddo
 927       enddo
 928       write (iout,'(/14a)') 'Residues',' nst ','  wst1',
 929      &  '    wst2','    wst3','    wst4','   dh11','   dh21',
 930      &  '   dh12','   dh22','    dt1','    dt2','    epsh1',
 931      &  '   sigh'
 932       do i=1,ntyp
 933         do j=1,i
 934           write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
 935      &      restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
 936      &       ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
 937      &       epshead(i,j),sig0head(i,j)
 938         enddo
 939       enddo
 940       write (iout,'(/12a)') 'Residues',' ch1',' ch2',
 941      &  '    rborn1','    rborn2','    wqdip1','    wqdip2',
 942      &  '     wquad'
 943       do i=1,ntyp
 944         do j=1,i
 945           write (iout,'(2(a3,1x),2i4,5f10.3)')
 946      &      restyp(i),restyp(j),icharge(i),icharge(j),
 947      &      rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
 948         enddo
 949       enddo
 950       write (iout,'(/12a)') 'Residues',
 951      &  '  alphpol1',
 952      &  '  alphpol2','  alphiso1','   alpiso2',
 953      &  '   alpiso3','   alpiso4','   sigiso1','   sigiso2',
 954      &  '     epsin'
 955       do i=1,ntyp
 956         do j=1,i
 957           write (iout,'(2(a3,1x),11f10.3)')
 958      &      restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
 959      &      (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
 960      &      epsintab(i,j)
 961         enddo
 962       enddo
 963       goto 70
 964
 965    60 continue
 966       close (isidep)
 967 C-----------------------------------------------------------------------
 968 C Calculate the "working" parameters of SC interactions.
 969
 970       IF (ipot.LT.6) THEN
 971        do i=1,ntyp
 972         do j=i,ntyp
 973          sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 974          sigma(j,i)=sigma(i,j)
 975          rs0(i,j)=dwa16*sigma(i,j)
 976          rs0(j,i)=rs0(i,j)
 977         enddo
 978        enddo
 979       END IF
 980
 981    70 continue
 982       write (iout,*) "IPOT=",ipot
 983       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 984      & 'Working parameters of the SC interactions:',
 985      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 986      & '  chi1   ','   chi2   '
 987       do i=1,ntyp
 988         do j=i,ntyp
 989           epsij=eps(i,j)
 990           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 ) THEN
 991             rrij=sigma(i,j)
 992           else
 993             rrij=rr0(i)+rr0(j)
 994           endif
 995           r0(i,j)=rrij
 996           r0(j,i)=rrij
 997           rrij=rrij**expon
 998           epsij=eps(i,j)
 999           sigeps=dsign(1.0D0,epsij)
1000           epsij=dabs(epsij)
1001           aa(i,j)=epsij*rrij*rrij
1002           bb(i,j)=-sigeps*epsij*rrij
1003           aa(j,i)=aa(i,j)
1004           bb(j,i)=bb(i,j)
1005           IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1006             sigt1sq=sigma0(i)**2
1007             sigt2sq=sigma0(j)**2
1008             sigii1=sigii(i)
1009             sigii2=sigii(j)
1010             ratsig1=sigt2sq/sigt1sq
1011             ratsig2=1.0D0/ratsig1
1012             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1013             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1014             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1015           else
1016             rsum_max=sigma(i,j)
1017           endif
1018 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1019             sigmaii(i,j)=rsum_max
1020             sigmaii(j,i)=rsum_max
1021 c         else
1022 c           sigmaii(i,j)=r0(i,j)
1023 c           sigmaii(j,i)=r0(i,j)
1024 c         endif
1025 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1026           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1027             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1028             augm(i,j)=epsij*r_augm**(2*expon)
1029 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1030             augm(j,i)=augm(i,j)
1031           else
1032             augm(i,j)=0.0D0
1033             augm(j,i)=0.0D0
1034           endif
1035           if (lprint) then
1036             if (ipot.lt.6) then
1037             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1038      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1039      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1040             else
1041             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1042      &       i3,40f10.4)')
1043      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1044      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1045      &      icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1046      &     (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1047      &     chis(i,j),chis(j,i),
1048      &     nstate(i,j),(wstate(k,i,j),k=1,4),
1049      &     ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1050      &     epshead(i,j),sig0head(i,j),
1051      &     rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1052      &     alphapol(i,j),alphapol(j,i),
1053      &     (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1054
1055             endif
1056                 endif
1057         enddo
1058       enddo
1059
1060 C
1061 C Define the SC-p interaction constants
1062 C
1063 #ifdef OLDSCP
1064       do i=1,20
1065 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1066 C helix formation)
1067 c       aad(i,1)=0.3D0*4.0D0**12
1068 C Following line for constants currently implemented
1069 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1070         aad(i,1)=1.5D0*4.0D0**12
1071 c       aad(i,1)=0.17D0*5.6D0**12
1072         aad(i,2)=aad(i,1)
1073 C "Soft" SC-p repulsion
1074         bad(i,1)=0.0D0
1075 C Following line for constants currently implemented
1076 c       aad(i,1)=0.3D0*4.0D0**6
1077 C "Hard" SC-p repulsion
1078         bad(i,1)=3.0D0*4.0D0**6
1079 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1080         bad(i,2)=bad(i,1)
1081 c       aad(i,1)=0.0D0
1082 c       aad(i,2)=0.0D0
1083 c       bad(i,1)=1228.8D0
1084 c       bad(i,2)=1228.8D0
1085       enddo
1086 #else
1087 C
1088 C 8/9/01 Read the SC-p interaction constants from file
1089 C
1090       do i=1,ntyp
1091         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1092       enddo
1093       do i=1,ntyp
1094         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1095         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1096         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1097         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1098       enddo
1099
1100       if (lprint) then
1101         write (iout,*) "Parameters of SC-p interactions:"
1102         do i=1,20
1103           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1104      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1105         enddo
1106       endif
1107 #endif
1108 C
1109 C Define the constants of the disulfide bridge
1110 C
1111       ebr=-5.50D0
1112 c
1113 c Old arbitrary potential - commented out.
1114 c
1115 c      dbr= 4.20D0
1116 c      fbr= 3.30D0
1117 c
1118 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1119 c energy surface of diethyl disulfide.
1120 c A. Liwo and U. Kozlowska, 11/24/03
1121 c
1122       D0CM = 3.78d0
1123       AKCM = 15.1d0
1124       AKTH = 11.0d0
1125       AKCT = 12.0d0
1126       V1SS =-1.08d0
1127       V2SS = 7.61d0
1128       V3SS = 13.7d0
1129
1130       if (lprint) then
1131       write (iout,'(/a)') "Disulfide bridge parameters:"
1132       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1133       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1134       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1135       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1136      & ' v3ss:',v3ss
1137       endif
1138       return
1139   110 write (iout,*) "Error reading bond energy parameters."
1140       goto 999
1141   111 write (iout,*) "Error reading bending energy parameters."
1142       goto 999
1143   112 write (iout,*) "Error reading rotamer energy parameters."
1144       goto 999
1145   113 write (iout,*) "Error reading torsional energy parameters."
1146       goto 999
1147   114 write (iout,*) "Error reading double torsional energy parameters."
1148       goto 999
1149   115 write (iout,*)
1150      &  "Error reading cumulant (multibody energy) parameters."
1151       goto 999
1152   116 write (iout,*) "Error reading electrostatic energy parameters."
1153       goto 999
1154   117 write (iout,*) "Error reading side chain interaction parameters."
1155       goto 999
1156   118 write (iout,*) "Error reading SCp interaction parameters."
1157       goto 999
1158   119 write (iout,*) "Error reading SCCOR parameters"
1159   999 continue
1160 #ifdef MPI
1161       call MPI_Finalize(Ierror)
1162 #endif
1163       stop
1164       end