source/wham/src-NEWSC/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,*) 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      &       l=1,maxinter)
 629             enddo
 630             write (iout,*) 'Lorenz constants'
 631             do k=1,nlor_sccor(i,j)
 632               write (iout,'(3(1pe15.5))')
 633      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 634             enddo
 635           enddo
 636         enddo
 637       endif
 638
 639 C
 640 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 641 C         interaction energy of the Gly, Ala, and Pro prototypes.
 642 C
 643       read (ifourier,*) nloctyp
 644       do i=1,nloctyp
 645         read (ifourier,*)
 646         read (ifourier,*) (b(ii,i),ii=1,13)
 647         if (lprint) then
 648         write (iout,*) 'Type',i
 649         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 650         endif
 651         B1(1,i)  = b(3,i)
 652         B1(2,i)  = b(5,i)
 653         B1tilde(1,i) = b(3,i)
 654         B1tilde(2,i) =-b(5,i)
 655         B2(1,i)  = b(2,i)
 656         B2(2,i)  = b(4,i)
 657         CC(1,1,i)= b(7,i)
 658         CC(2,2,i)=-b(7,i)
 659         CC(2,1,i)= b(9,i)
 660         CC(1,2,i)= b(9,i)
 661         Ctilde(1,1,i)=b(7,i)
 662         Ctilde(1,2,i)=b(9,i)
 663         Ctilde(2,1,i)=-b(9,i)
 664         Ctilde(2,2,i)=b(7,i)
 665         DD(1,1,i)= b(6,i)
 666         DD(2,2,i)=-b(6,i)
 667         DD(2,1,i)= b(8,i)
 668         DD(1,2,i)= b(8,i)
 669         Dtilde(1,1,i)=b(6,i)
 670         Dtilde(1,2,i)=b(8,i)
 671         Dtilde(2,1,i)=-b(8,i)
 672         Dtilde(2,2,i)=b(6,i)
 673         EE(1,1,i)= b(10,i)+b(11,i)
 674         EE(2,2,i)=-b(10,i)+b(11,i)
 675         EE(2,1,i)= b(12,i)-b(13,i)
 676         EE(1,2,i)= b(12,i)+b(13,i)
 677       enddo
 678       if (lprint) then
 679       do i=1,nloctyp
 680         write (iout,*) 'Type',i
 681         write (iout,*) 'B1'
 682 c        write (iout,'(f10.5)') B1(:,i)
 683         write(iout,*) B1(1,i),B1(2,i)
 684         write (iout,*) 'B2'
 685 c        write (iout,'(f10.5)') B2(:,i)
 686         write(iout,*) B2(1,i),B2(2,i)
 687         write (iout,*) 'CC'
 688         do j=1,2
 689           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 690         enddo
 691         write(iout,*) 'DD'
 692         do j=1,2
 693           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 694         enddo
 695         write(iout,*) 'EE'
 696         do j=1,2
 697           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 698         enddo
 699       enddo
 700       endif
 701 C
 702 C Read electrostatic-interaction parameters
 703 C
 704       if (lprint) then
 705         write (iout,'(/a)') 'Electrostatic interaction constants:'
 706         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 707      &            'IT','JT','APP','BPP','AEL6','AEL3'
 708       endif
 709       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 710       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 711       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 712       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 713       close (ielep)
 714       do i=1,2
 715         do j=1,2
 716         rri=rpp(i,j)**6
 717         app (i,j)=epp(i,j)*rri*rri
 718         bpp (i,j)=-2.0D0*epp(i,j)*rri
 719         ael6(i,j)=elpp6(i,j)*4.2D0**6
 720         ael3(i,j)=elpp3(i,j)*4.2D0**3
 721         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 722      &                    ael6(i,j),ael3(i,j)
 723         enddo
 724       enddo
 725 C
 726 C Read side-chain interaction parameters.
 727 C
 728       read (isidep,*) ipot,expon
 729       if (ipot.lt.1 .or. ipot.gt.6) then
 730         write (iout,'(2a)') 'Error while reading SC interaction',
 731      &               'potential file - unknown potential type.'
 732         stop
 733       endif
 734       expon2=expon/2
 735       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 736      & ', exponents are ',expon,2*expon
 737       goto (10,20,30,30,40,50) ipot
 738 C----------------------- LJ potential ---------------------------------
 739    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 740       if (lprint) then
 741         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 742         write (iout,'(a/)') 'The epsilon array:'
 743         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 744         write (iout,'(/a)') 'One-body parameters:'
 745         write (iout,'(a,4x,a)') 'residue','sigma'
 746         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 747       endif
 748       goto 60
 749 C----------------------- LJK potential --------------------------------
 750    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 751      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 752       if (lprint) then
 753         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 754         write (iout,'(a/)') 'The epsilon array:'
 755         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 756         write (iout,'(/a)') 'One-body parameters:'
 757         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 758         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 759      &        i=1,ntyp)
 760       endif
 761       goto 60
 762 C---------------------- GB or BP potential -----------------------------
 763    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 764      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 765      &  (alp(i),i=1,ntyp)
 766 C For the GB potential convert sigma'**2 into chi'
 767       if (ipot.eq.4) then
 768         do i=1,ntyp
 769           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 770         enddo
 771       endif
 772       if (lprint) then
 773         write (iout,'(/a/)') 'Parameters of the BP potential:'
 774         write (iout,'(a/)') 'The epsilon array:'
 775         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 776         write (iout,'(/a)') 'One-body parameters:'
 777         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 778      &       '    chip  ','    alph  '
 779         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 780      &                     chip(i),alp(i),i=1,ntyp)
 781       endif
 782       goto 60
 783 C--------------------- GBV potential -----------------------------------
 784    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 785      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 786      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 787       if (lprint) then
 788         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 789         write (iout,'(a/)') 'The epsilon array:'
 790         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 791         write (iout,'(/a)') 'One-body parameters:'
 792         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 793      &      's||/s_|_^2','    chip  ','    alph  '
 794         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 795      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 796       endif
 797       goto 60
 798 C--------------------- Momo potential -----------------------------------
 799
 800    50 continue
 801
 802       read (isidep,*) (icharge(i),i=1,ntyp)
 803 c      write (2,*) "icharge",(icharge(i),i=1,ntyp)
 804       do i=1,ntyp
 805        do j=1,i
 806 c!        write (*,*) "Im in ", i, " ", j
 807         read(isidep,*)
 808      &  eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
 809      &  (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
 810      &  chis(i,j),chis(j,i),
 811      &  nstate(i,j),(wstate(k,i,j),k=1,4),
 812      &  dhead(1,1,i,j),
 813      &  dhead(1,2,i,j),
 814      &  dhead(2,1,i,j),
 815      &  dhead(2,2,i,j),
 816      &  dtail(1,i,j),dtail(2,i,j),
 817      &  epshead(i,j),sig0head(i,j),
 818      &  rborn(i,j),rborn(j,i),
 819      &  (wqdip(k,i,j),k=1,2),wquad(i,j),
 820      &  alphapol(i,j),alphapol(j,i),
 821      &  (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
 822        END DO
 823       END DO
 824 c!      write (*,*) "nstate gly-gly", nstate(10,10)
 825 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
 826 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
 827 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
 828 c! FROM HIS-ARG.
 829 c!
 830 c! DISABLE IT AT >>YOUR OWN PERIL<<
 831 c!
 832       DO i = 1, ntyp
 833        DO j = i+1, ntyp
 834         eps(i,j) = eps(j,i)
 835         sigma(i,j) = sigma(j,i)
 836         nstate(i,j) = nstate(j,i)
 837         sigmap1(i,j) = sigmap1(j,i)
 838         sigmap2(i,j) = sigmap2(j,i)
 839         sigiso1(i,j) = sigiso1(j,i)
 840         sigiso2(i,j) = sigiso2(j,i)
 841
 842         DO k = 1, 4
 843          alphasur(k,i,j) = alphasur(k,j,i)
 844          wstate(k,i,j) = wstate(k,j,i)
 845          alphiso(k,i,j) = alphiso(k,j,i)
 846         END DO
 847
 848         dhead(2,1,i,j) = dhead(1,1,j,i)
 849         dhead(2,2,i,j) = dhead(1,2,j,i)
 850         dhead(1,1,i,j) = dhead(2,1,j,i)
 851         dhead(1,2,i,j) = dhead(2,2,j,i)
 852         dtail(1,i,j) = dtail(1,j,i)
 853         dtail(2,i,j) = dtail(2,j,i)
 854 c!        DO k = 1, 2
 855 c!         DO l = 1, 2
 856 c!         write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
 857 c!         write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
 858 c!          dhead(l,k,i,j) = dhead(k,l,j,i)
 859 c!         END DO
 860 c!        END DO
 861
 862         epshead(i,j) = epshead(j,i)
 863         sig0head(i,j) = sig0head(j,i)
 864
 865         DO k = 1, 2
 866          wqdip(k,i,j) = wqdip(k,j,i)
 867         END DO
 868
 869         wquad(i,j) = wquad(j,i)
 870         epsintab(i,j) = epsintab(j,i)
 871
 872        END DO
 873       END DO
 874
 875       if (.not.lprint) goto 70
 876       write (iout,'(a)')
 877      & "Parameters of the new physics-based SC-SC interaction potential"
 878       write (iout,'(/7a)') 'Residues','     epsGB','       rGB',
 879      &  '    chi1GB','    chi2GB','   chip1GB','   chip2GB'
 880       do i=1,ntyp
 881         do j=1,i
 882           write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
 883      &      restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
 884      &       chipp(i,j),chipp(j,i)
 885         enddo
 886       enddo
 887       write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
 888      &  ' alphasur3',' alphasur4','   sigmap1','   sigmap2',
 889      &  '     chis1','     chis2'
 890       do i=1,ntyp
 891         do j=1,i
 892           write (iout,'(2(a3,1x),8(0pf10.3))')
 893      &      restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
 894      &       sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
 895         enddo
 896       enddo
 897       write (iout,'(/14a)') 'Residues',' nst ','  wst1',
 898      &  '    wst2','    wst3','    wst4','   dh11','   dh21',
 899      &  '   dh12','   dh22','    dt1','    dt2','    epsh1',
 900      &  '   sigh'
 901       do i=1,ntyp
 902         do j=1,i
 903           write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
 904      &      restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
 905      &       ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
 906      &       epshead(i,j),sig0head(i,j)
 907         enddo
 908       enddo
 909       write (iout,'(/12a)') 'Residues',' ch1',' ch2',
 910      &  '    rborn1','    rborn2','    wqdip1','    wqdip2',
 911      &  '     wquad'
 912       do i=1,ntyp
 913         do j=1,i
 914           write (iout,'(2(a3,1x),2i4,5f10.3)')
 915      &      restyp(i),restyp(j),icharge(i),icharge(j),
 916      &      rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
 917         enddo
 918       enddo
 919       write (iout,'(/12a)') 'Residues',
 920      &  '  alphpol1',
 921      &  '  alphpol2','  alphiso1','   alpiso2',
 922      &  '   alpiso3','   alpiso4','   sigiso1','   sigiso2',
 923      &  '     epsin'
 924       do i=1,ntyp
 925         do j=1,i
 926           write (iout,'(2(a3,1x),11f10.3)')
 927      &      restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
 928      &      (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
 929      &      epsintab(i,j)
 930         enddo
 931       enddo
 932       goto 70
 933
 934    60 continue
 935       close (isidep)
 936 C-----------------------------------------------------------------------
 937 C Calculate the "working" parameters of SC interactions.
 938
 939       IF (ipot.LT.6) THEN
 940        do i=1,ntyp
 941         do j=i,ntyp
 942          sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 943          sigma(j,i)=sigma(i,j)
 944          rs0(i,j)=dwa16*sigma(i,j)
 945          rs0(j,i)=rs0(i,j)
 946         enddo
 947        enddo
 948       END IF
 949
 950    70 continue
 951       write (iout,*) "IPOT=",ipot
 952       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 953      & 'Working parameters of the SC interactions:',
 954      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 955      & '  chi1   ','   chi2   '
 956       do i=1,ntyp
 957         do j=i,ntyp
 958           epsij=eps(i,j)
 959           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 ) THEN
 960             rrij=sigma(i,j)
 961           else
 962             rrij=rr0(i)+rr0(j)
 963           endif
 964           r0(i,j)=rrij
 965           r0(j,i)=rrij
 966           rrij=rrij**expon
 967           epsij=eps(i,j)
 968           sigeps=dsign(1.0D0,epsij)
 969           epsij=dabs(epsij)
 970           aa(i,j)=epsij*rrij*rrij
 971           bb(i,j)=-sigeps*epsij*rrij
 972           aa(j,i)=aa(i,j)
 973           bb(j,i)=bb(i,j)
 974           IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
 975             sigt1sq=sigma0(i)**2
 976             sigt2sq=sigma0(j)**2
 977             sigii1=sigii(i)
 978             sigii2=sigii(j)
 979             ratsig1=sigt2sq/sigt1sq
 980             ratsig2=1.0D0/ratsig1
 981             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 982             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 983             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 984           else
 985             rsum_max=sigma(i,j)
 986           endif
 987 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 988             sigmaii(i,j)=rsum_max
 989             sigmaii(j,i)=rsum_max
 990 c         else
 991 c           sigmaii(i,j)=r0(i,j)
 992 c           sigmaii(j,i)=r0(i,j)
 993 c         endif
 994 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 995           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 996             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 997             augm(i,j)=epsij*r_augm**(2*expon)
 998 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 999             augm(j,i)=augm(i,j)
1000           else
1001             augm(i,j)=0.0D0
1002             augm(j,i)=0.0D0
1003           endif
1004           if (lprint) then
1005             if (ipot.lt.6) then
1006             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1007      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1008      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1009             else
1010             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1011      &       i3,40f10.4)')
1012      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1013      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1014      &      icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1015      &     (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1016      &     chis(i,j),chis(j,i),
1017      &     nstate(i,j),(wstate(k,i,j),k=1,4),
1018      &     ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1019      &     epshead(i,j),sig0head(i,j),
1020      &     rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1021      &     alphapol(i,j),alphapol(j,i),
1022      &     (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1023
1024             endif
1025                 endif
1026         enddo
1027       enddo
1028
1029 C
1030 C Define the SC-p interaction constants
1031 C
1032 #ifdef OLDSCP
1033       do i=1,20
1034 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1035 C helix formation)
1036 c       aad(i,1)=0.3D0*4.0D0**12
1037 C Following line for constants currently implemented
1038 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1039         aad(i,1)=1.5D0*4.0D0**12
1040 c       aad(i,1)=0.17D0*5.6D0**12
1041         aad(i,2)=aad(i,1)
1042 C "Soft" SC-p repulsion
1043         bad(i,1)=0.0D0
1044 C Following line for constants currently implemented
1045 c       aad(i,1)=0.3D0*4.0D0**6
1046 C "Hard" SC-p repulsion
1047         bad(i,1)=3.0D0*4.0D0**6
1048 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1049         bad(i,2)=bad(i,1)
1050 c       aad(i,1)=0.0D0
1051 c       aad(i,2)=0.0D0
1052 c       bad(i,1)=1228.8D0
1053 c       bad(i,2)=1228.8D0
1054       enddo
1055 #else
1056 C
1057 C 8/9/01 Read the SC-p interaction constants from file
1058 C
1059       do i=1,ntyp
1060         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1061       enddo
1062       do i=1,ntyp
1063         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1064         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1065         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1066         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1067       enddo
1068
1069       if (lprint) then
1070         write (iout,*) "Parameters of SC-p interactions:"
1071         do i=1,20
1072           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1073      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1074         enddo
1075       endif
1076 #endif
1077 C
1078 C Define the constants of the disulfide bridge
1079 C
1080       ebr=-5.50D0
1081 c
1082 c Old arbitrary potential - commented out.
1083 c
1084 c      dbr= 4.20D0
1085 c      fbr= 3.30D0
1086 c
1087 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1088 c energy surface of diethyl disulfide.
1089 c A. Liwo and U. Kozlowska, 11/24/03
1090 c
1091       D0CM = 3.78d0
1092       AKCM = 15.1d0
1093       AKTH = 11.0d0
1094       AKCT = 12.0d0
1095       V1SS =-1.08d0
1096       V2SS = 7.61d0
1097       V3SS = 13.7d0
1098
1099       if (lprint) then
1100       write (iout,'(/a)') "Disulfide bridge parameters:"
1101       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1102       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1103       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1104       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1105      & ' v3ss:',v3ss
1106       endif
1107       return
1108       end