source/wham/src-M-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       character*1 toronelet(-2:2) /"p","a","G","A","P"/
  27       logical lprint
  28       dimension blower(3,3,maxlob)
  29       character*800 controlcard
  30       character*256 bondname_t,thetname_t,rotname_t,torname_t,
  31      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
  32      & sccorname_t
  33       integer ilen
  34       external ilen
  35       character*16 key
  36       integer iparm
  37       double precision ip,mp
  38 C
  39 C Body
  40 C
  41 C Set LPRINT=.TRUE. for debugging
  42       dwa16=2.0d0**(1.0d0/6.0d0)
  43       lprint=.true.
  44       itypro=20
  45 C Assign virtual-bond length
  46       vbl=3.8D0
  47       vblinv=1.0D0/vbl
  48       vblinv2=vblinv*vblinv
  49       call card_concat(controlcard,.true.)
  50       wname(4)="WCORRH"
  51       do i=1,n_ene
  52         key = wname(i)(:ilen(wname(i)))
  53         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
  54       enddo
  55
  56       write (iout,*) "iparm",iparm," myparm",myparm
  57 c If reading not own parameters, skip assignment
  58
  59       if (iparm.eq.myparm .or. .not.separate_parset) then
  60
  61 c
  62 c Setup weights for UNRES
  63 c
  64       wsc=ww(1)
  65       wscp=ww(2)
  66       welec=ww(3)
  67       wcorr=ww(4)
  68       wcorr5=ww(5)
  69       wcorr6=ww(6)
  70       wel_loc=ww(7)
  71       wturn3=ww(8)
  72       wturn4=ww(9)
  73       wturn6=ww(10)
  74       wang=ww(11)
  75       wscloc=ww(12)
  76       wtor=ww(13)
  77       wtor_d=ww(14)
  78       wvdwpp=ww(16)
  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,"SCCORPAR",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,1,1),j=1,2),
 198      &    (bthet(j,i,1,1),j=1,2)
 199         read (ithep,*) (polthet(j,i),j=0,3)
 200         read (ithep,*) (gthet(j,i),j=1,3)
 201         read (ithep,*) theta0(i),sig0(i),sigc0(i)
 202         sigc0(i)=sigc0(i)**2
 203       enddo
 204       do i=1,ntyp
 205       athet(1,i,1,-1)=athet(1,i,1,1)
 206       athet(2,i,1,-1)=athet(2,i,1,1)
 207       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 208       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 209       athet(1,i,-1,1)=-athet(1,i,1,1)
 210       athet(2,i,-1,1)=-athet(2,i,1,1)
 211       bthet(1,i,-1,1)=bthet(1,i,1,1)
 212       bthet(2,i,-1,1)=bthet(2,i,1,1)
 213       enddo
 214       do i=-ntyp,-1
 215       a0thet(i)=a0thet(-i)
 216       athet(1,i,-1,-1)=athet(1,-i,1,1)
 217       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 218       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 219       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 220       athet(1,i,-1,1)=athet(1,-i,1,1)
 221       athet(2,i,-1,1)=-athet(2,-i,1,1)
 222       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 223       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 224       athet(1,i,1,-1)=-athet(1,-i,1,1)
 225       athet(2,i,1,-1)=athet(2,-i,1,1)
 226       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 227       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 228       theta0(i)=theta0(-i)
 229       sig0(i)=sig0(-i)
 230       sigc0(i)=sigc0(-i)
 231        do j=0,3
 232         polthet(j,i)=polthet(j,-i)
 233        enddo
 234        do j=1,3
 235          gthet(j,i)=gthet(j,-i)
 236        enddo
 237       enddo
 238       close (ithep)
 239       if (lprint) then
 240 c       write (iout,'(a)')
 241 c    &    'Parameters of the virtual-bond valence angles:'
 242 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 243 c    & '    ATHETA0   ','         A1   ','        A2    ',
 244 c    & '        B1    ','         B2   '
 245 c       do i=1,ntyp
 246 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 247 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 248 c       enddo
 249 c       write (iout,'(/a/9x,5a/79(1h-))')
 250 c    & 'Parameters of the expression for sigma(theta_c):',
 251 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 252 c    & '     ALPH3    ','    SIGMA0C   '
 253 c       do i=1,ntyp
 254 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 255 c    &      (polthet(j,i),j=0,3),sigc0(i)
 256 c       enddo
 257 c       write (iout,'(/a/9x,5a/79(1h-))')
 258 c    & 'Parameters of the second gaussian:',
 259 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 260 c    & '        G2    ','         G3   '
 261 c       do i=1,ntyp
 262 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 263 c    &       sig0(i),(gthet(j,i),j=1,3)
 264 c       enddo
 265         write (iout,'(a)')
 266      &    'Parameters of the virtual-bond valence angles:'
 267         write (iout,'(/a/9x,5a/79(1h-))')
 268      & 'Coefficients of expansion',
 269      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 270      & '   b1*10^1    ','    b2*10^1   '
 271         do i=1,ntyp
 272           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 273      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
 274      &        (10*bthet(j,i,1,1),j=1,2)
 275         enddo
 276         write (iout,'(/a/9x,5a/79(1h-))')
 277      & 'Parameters of the expression for sigma(theta_c):',
 278      & ' alpha0       ','  alph1       ',' alph2        ',
 279      & ' alhp3        ','   sigma0c    '
 280         do i=1,ntyp
 281           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 282      &      (polthet(j,i),j=0,3),sigc0(i)
 283         enddo
 284         write (iout,'(/a/9x,5a/79(1h-))')
 285      & 'Parameters of the second gaussian:',
 286      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 287      & '        G2    ','   G3*10^1    '
 288         do i=1,ntyp
 289           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 290      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 291         enddo
 292       endif
 293 #else
 294 C
 295 C Read the parameters of Utheta determined from ab initio surfaces
 296 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 297 C
 298 c      write (iout,*) "tu dochodze"
 299       read (ithep,*) nthetyp,ntheterm,ntheterm2,
 300      &  ntheterm3,nsingle,ndouble
 301       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 302       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 303       do i=-ntyp1,-1
 304         ithetyp(i)=-ithetyp(-i)
 305       enddo
 306 c      write (iout,*) "tu dochodze"
 307       do iblock=1,2
 308       do i=-maxthetyp,maxthetyp
 309         do j=-maxthetyp,maxthetyp
 310           do k=-maxthetyp,maxthetyp
 311             aa0thet(i,j,k,iblock)=0.0d0
 312             do l=1,ntheterm
 313               aathet(l,i,j,k,iblock)=0.0d0
 314             enddo
 315             do l=1,ntheterm2
 316               do m=1,nsingle
 317                 bbthet(m,l,i,j,k,iblock)=0.0d0
 318                 ccthet(m,l,i,j,k,iblock)=0.0d0
 319                 ddthet(m,l,i,j,k,iblock)=0.0d0
 320                 eethet(m,l,i,j,k,iblock)=0.0d0
 321               enddo
 322             enddo
 323             do l=1,ntheterm3
 324               do m=1,ndouble
 325                 do mm=1,ndouble
 326                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 327                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 328                 enddo
 329               enddo
 330             enddo
 331           enddo
 332         enddo
 333       enddo
 334       enddo
 335       do iblock=1,2
 336       do i=0,nthetyp
 337         do j=-nthetyp,nthetyp
 338           do k=-nthetyp,nthetyp
 339             read (ithep,'(6a)') res1
 340             read (ithep,*) aa0thet(i,j,k,iblock)
 341             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 342             read (ithep,*)
 343      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 344      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 345      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 346      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 347      &        ,ll=1,ntheterm2)
 348             read (ithep,*)
 349      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 350      &      ffthet(lll,llll,ll,i,j,k,iblock),
 351      &         ggthet(llll,lll,ll,i,j,k,iblock)
 352      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 353      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 354           enddo
 355         enddo
 356       enddo
 357 C
 358 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 359 C coefficients of theta-and-gamma-dependent terms are zero.
 360 C
 361       do i=1,nthetyp
 362         do j=1,nthetyp
 363           do l=1,ntheterm
 364             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 365             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 366           enddo
 367           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 368           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 369         enddo
 370         do l=1,ntheterm
 371           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 372         enddo
 373         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 374       enddo
 375       enddo
 376 C Substitution for D aminoacids from symmetry.
 377       do iblock=1,2
 378       do i=-nthetyp,0
 379         do j=-nthetyp,nthetyp
 380           do k=-nthetyp,nthetyp
 381            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 382            do l=1,ntheterm
 383            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 384            enddo
 385            do ll=1,ntheterm2
 386             do lll=1,nsingle
 387             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 388             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 389             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 390             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 391             enddo
 392           enddo
 393           do ll=1,ntheterm3
 394            do lll=2,ndouble
 395             do llll=1,lll-1
 396             ffthet(llll,lll,ll,i,j,k,iblock)=
 397      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 398             ffthet(lll,llll,ll,i,j,k,iblock)=
 399      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 400             ggthet(llll,lll,ll,i,j,k,iblock)=
 401      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 402             ggthet(lll,llll,ll,i,j,k,iblock)=
 403      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 404             enddo !ll
 405            enddo  !lll
 406           enddo   !llll
 407          enddo    !k
 408         enddo     !j
 409        enddo      !i
 410       enddo       !iblock
 411
 412 C
 413 C Control printout of the coefficients of virtual-bond-angle potentials
 414 C
 415       if (lprint) then
 416         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 417         do iblock=1,2
 418         do i=1,nthetyp+1
 419           do j=1,nthetyp+1
 420             do k=1,nthetyp+1
 421               write (iout,'(//4a)')
 422      &         'Type ',onelett(i),onelett(j),onelett(k)
 423               write (iout,'(//a,10x,a)') " l","a[l]"
 424               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 425               write (iout,'(i2,1pe15.5)')
 426      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 427             do l=1,ntheterm2
 428               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 429      &          "b",l,"c",l,"d",l,"e",l
 430               do m=1,nsingle
 431                 write (iout,'(i2,4(1pe15.5))') m,
 432      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 433      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 434               enddo
 435             enddo
 436             do l=1,ntheterm3
 437               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 438      &          "f+",l,"f-",l,"g+",l,"g-",l
 439               do m=2,ndouble
 440                 do n=1,m-1
 441                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 442      &              ffthet(n,m,l,i,j,k,iblock),
 443      &              ffthet(m,n,l,i,j,k,iblock),
 444      &              ggthet(n,m,l,i,j,k,iblock),
 445      &              ggthet(m,n,l,i,j,k,iblock)
 446                 enddo
 447               enddo
 448             enddo
 449           enddo
 450         enddo
 451         enddo
 452       enddo
 453       call flush(iout)
 454       endif
 455 #endif
 456
 457 #ifdef CRYST_SC
 458 C
 459 C Read the parameters of the probability distribution/energy expression
 460 C of the side chains.
 461 C
 462       do i=1,ntyp
 463 cc      write (iout,*) "tu dochodze",i
 464         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 465         if (i.eq.10) then
 466           dsc_inv(i)=0.0D0
 467         else
 468           dsc_inv(i)=1.0D0/dsc(i)
 469         endif
 470         if (i.ne.10) then
 471         do j=1,nlob(i)
 472           do k=1,3
 473             do l=1,3
 474               blower(l,k,j)=0.0D0
 475             enddo
 476           enddo
 477         enddo
 478         bsc(1,i)=0.0D0
 479         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 480         censc(1,1,-i)=censc(1,1,i)
 481         censc(2,1,-i)=censc(2,1,i)
 482         censc(3,1,-i)=-censc(3,1,i)
 483         do j=2,nlob(i)
 484           read (irotam,*) bsc(j,i)
 485           read (irotam,*) (censc(k,j,i),k=1,3),
 486      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 487         censc(1,j,-i)=censc(1,j,i)
 488         censc(2,j,-i)=censc(2,j,i)
 489         censc(3,j,-i)=-censc(3,j,i)
 490 C BSC is amplitude of Gaussian
 491         enddo
 492         do j=1,nlob(i)
 493           do k=1,3
 494             do l=1,k
 495               akl=0.0D0
 496               do m=1,3
 497                 akl=akl+blower(k,m,j)*blower(l,m,j)
 498               enddo
 499               gaussc(k,l,j,i)=akl
 500               gaussc(l,k,j,i)=akl
 501              if (((k.eq.3).and.(l.ne.3))
 502      &        .or.((l.eq.3).and.(k.ne.3))) then
 503                 gaussc(k,l,j,-i)=-akl
 504                 gaussc(l,k,j,-i)=-akl
 505               else
 506                 gaussc(k,l,j,-i)=akl
 507                 gaussc(l,k,j,-i)=akl
 508               endif
 509             enddo
 510           enddo
 511         enddo
 512         endif
 513       enddo
 514       close (irotam)
 515       if (lprint) then
 516         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 517         do i=1,ntyp
 518           nlobi=nlob(i)
 519           if (nlobi.gt.0) then
 520           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 521      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 522 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 523 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 524 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 525 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 526            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 527      &                             'log h',(bsc(j,i),j=1,nlobi)
 528            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 529      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 530 c          write (iout,'(a)')
 531 c         do j=1,nlobi
 532 c           ind=0
 533 c           do k=1,3
 534 c             do l=1,k
 535 c              ind=ind+1
 536 c              blower(k,l,j)=gaussc(ind,j,i)
 537 c             enddo
 538 c           enddo
 539 c         enddo
 540           do k=1,3
 541             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 542      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 543           enddo
 544           endif
 545         enddo
 546       endif
 547 #else
 548 C
 549 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 550 C added by Urszula Kozlowska 07/11/2007
 551 C
 552       do i=1,ntyp
 553         read (irotam,*)
 554        if (i.eq.10) then
 555          read (irotam,*)
 556        else
 557          do j=1,65
 558            read(irotam,*) sc_parmin(j,i)
 559          enddo
 560        endif
 561       enddo
 562 #endif
 563       close(irotam)
 564 #ifdef CRYST_TOR
 565 C
 566 C Read torsional parameters in old format
 567 C
 568       read (itorp,*) ntortyp,nterm_old
 569       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 570       read (itorp,*) (itortyp(i),i=1,ntyp)
 571       do i=1,ntortyp
 572         do j=1,ntortyp
 573           read (itorp,'(a)')
 574           do k=1,nterm_old
 575             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 576           enddo
 577         enddo
 578       enddo
 579       close (itorp)
 580       if (lprint) then
 581         write (iout,'(/a/)') 'Torsional constants:'
 582         do i=1,ntortyp
 583           do j=1,ntortyp
 584             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 585             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 586           enddo
 587         enddo
 588       endif
 589
 590
 591 #else
 592 C
 593 C Read torsional parameters
 594 C
 595       read (itorp,*) ntortyp
 596       read (itorp,*) (itortyp(i),i=1,ntyp)
 597       write (iout,*) 'ntortyp',ntortyp
 598       do iblock=1,2
 599       do i=-ntyp,-1
 600        itortyp(i)=-itortyp(-i)
 601       enddo
 602 c      write (iout,*) 'ntortyp',ntortyp
 603       do i=0,ntortyp-1
 604         do j=-ntortyp+1,ntortyp-1
 605           read (itorp,*) nterm(i,j,iblock),
 606      &          nlor(i,j,iblock)
 607           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 608           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 609           v0ij=0.0d0
 610           si=-1.0d0
 611           do k=1,nterm(i,j,iblock)
 612             read (itorp,*) kk,v1(k,i,j,iblock),
 613      &      v2(k,i,j,iblock)
 614             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 615             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 616             v0ij=v0ij+si*v1(k,i,j,iblock)
 617             si=-si
 618          enddo
 619           do k=1,nlor(i,j,iblock)
 620             read (itorp,*) kk,vlor1(k,i,j),
 621      &        vlor2(k,i,j),vlor3(k,i,j)
 622             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 623           enddo
 624           v0(i,j,iblock)=v0ij
 625           v0(-i,-j,iblock)=v0ij
 626         enddo
 627       enddo
 628       enddo
 629       close (itorp)
 630       if (lprint) then
 631         write (iout,'(/a/)') 'Torsional constants:'
 632         do i=1,ntortyp
 633           do j=1,ntortyp
 634             write (iout,*) 'ityp',i,' jtyp',j
 635             write (iout,*) 'Fourier constants'
 636             do iblock=1,2
 637             do k=1,nterm(i,j,iblock)
 638               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 639      &        v2(k,i,j,iblock)
 640             enddo
 641             write (iout,*) 'Lorenz constants'
 642             do k=1,nlor(i,j,iblock)
 643               write (iout,'(3(1pe15.5))')
 644      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 645             enddo
 646             enddo
 647           enddo
 648         enddo
 649       endif
 650 C
 651 C 6/23/01 Read parameters for double torsionals
 652 C
 653       do iblock=1,2
 654       do i=0,ntortyp-1
 655         do j=-ntortyp+1,ntortyp-1
 656           do k=-ntortyp+1,ntortyp-1
 657             read (itordp,'(3a1)') t1,t2,t3
 658 c              write (iout,*) "OK onelett",
 659 c     &         i,j,k,t1,t2,t3
 660
 661             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 662      &        .or. t3.ne.toronelet(k)) then
 663               write (iout,*) "Error in double torsional parameter file",
 664      &         i,j,k,t1,t2,t3
 665 #ifdef MPI
 666               call MPI_Finalize(Ierror)
 667 #endif
 668                stop "Error in double torsional parameter file"
 669             endif
 670           read (itordp,*) ntermd_1(i,j,k,iblock),
 671      &         ntermd_2(i,j,k,iblock)
 672             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 673             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 674             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
 675      &         ntermd_1(i,j,k,iblock))
 676             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
 677      &         ntermd_1(i,j,k,iblock))
 678             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
 679      &         ntermd_1(i,j,k,iblock))
 680             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
 681      &         ntermd_1(i,j,k,iblock))
 682 C Martix of D parameters for one dimesional foureir series
 683             do l=1,ntermd_1(i,j,k,iblock)
 684              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 685              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 686              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 687              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 688 c            write(iout,*) "whcodze" ,
 689 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 690             enddo
 691             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
 692      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 693      &         v2s(m,l,i,j,k,iblock),
 694      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 695 C Martix of D parameters for two dimesional fourier series
 696             do l=1,ntermd_2(i,j,k,iblock)
 697              do m=1,l-1
 698              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 699              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 700              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 701              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 702              enddo!m
 703             enddo!l
 704           enddo!k
 705         enddo!j
 706       enddo!i
 707       enddo!iblock
 708       if (lprint) then
 709       write (iout,*)
 710       write (iout,*) 'Constants for double torsionals'
 711       do iblock=1,2
 712       do i=0,ntortyp-1
 713         do j=-ntortyp+1,ntortyp-1
 714           do k=-ntortyp+1,ntortyp-1
 715             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 716      &        ' nsingle',ntermd_1(i,j,k,iblock),
 717      &        ' ndouble',ntermd_2(i,j,k,iblock)
 718             write (iout,*)
 719             write (iout,*) 'Single angles:'
 720             do l=1,ntermd_1(i,j,k,iblock)
 721               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 722      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 723      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 724      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 725             enddo
 726             write (iout,*)
 727             write (iout,*) 'Pairs of angles:'
 728             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 729             do l=1,ntermd_2(i,j,k,iblock)
 730               write (iout,'(i5,20f10.5)')
 731      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 732             enddo
 733             write (iout,*)
 734            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 735             do l=1,ntermd_2(i,j,k,iblock)
 736               write (iout,'(i5,20f10.5)')
 737      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 738      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 739             enddo
 740             write (iout,*)
 741           enddo
 742         enddo
 743       enddo
 744       enddo
 745       endif
 746 #endif
 747 C Read of Side-chain backbone correlation parameters
 748 C Modified 11 May 2012 by Adasko
 749 CCC
 750 C
 751       read (isccor,*) nsccortyp
 752       read (isccor,*) (isccortyp(i),i=1,ntyp)
 753       do i=-ntyp,-1
 754         isccortyp(i)=-isccortyp(-i)
 755       enddo
 756       iscprol=isccortyp(20)
 757 c      write (iout,*) 'ntortyp',ntortyp
 758       maxinter=3
 759 cc maxinter is maximum interaction sites
 760       do l=1,maxinter
 761       do i=1,nsccortyp
 762         do j=1,nsccortyp
 763           read (isccor,*)
 764      &nterm_sccor(i,j),nlor_sccor(i,j)
 765           write (iout,*) nterm_sccor(i,j)
 766           v0ijsccor=0.0d0
 767           v0ijsccor1=0.0d0
 768           v0ijsccor2=0.0d0
 769           v0ijsccor3=0.0d0
 770           si=-1.0d0
 771           nterm_sccor(-i,j)=nterm_sccor(i,j)
 772           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 773           nterm_sccor(i,-j)=nterm_sccor(i,j)
 774           write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
 775      &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
 776           do k=1,nterm_sccor(i,j)
 777             read (isccor,*) kk,v1sccor(k,l,i,j)
 778      &    ,v2sccor(k,l,i,j)
 779             if (j.eq.iscprol) then
 780              if (i.eq.isccortyp(10)) then
 781              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 782              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 783              else
 784              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 785      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 786              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 787      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 788              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 789              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 790              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 791              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 792              endif
 793             else
 794              if (i.eq.isccortyp(10)) then
 795              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 796              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 797              else
 798                if (j.eq.isccortyp(10)) then
 799              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 800              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 801                else
 802              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 803              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 804              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 805              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 806              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 807              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 808                 endif
 809                endif
 810             endif
 811             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 812             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 813             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 814             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 815             si=-si
 816            enddo
 817           do k=1,nlor_sccor(i,j)
 818             read (isccor,*) kk,vlor1sccor(k,i,j),
 819      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 820             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 821      &(1+vlor3sccor(k,i,j)**2)
 822           enddo
 823           v0sccor(l,i,j)=v0ijsccor
 824           v0sccor(l,-i,j)=v0ijsccor1
 825           v0sccor(l,i,-j)=v0ijsccor2
 826           v0sccor(l,-i,-j)=v0ijsccor3
 827           enddo
 828         enddo
 829       enddo
 830       close (isccor)
 831       if (lprint) then
 832         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 833         do i=1,nsccortyp
 834           do j=1,nsccortyp
 835             write (iout,*) 'ityp',i,' jtyp',j
 836             write (iout,*) 'Fourier constants'
 837             do l=1,3
 838             do k=1,nterm_sccor(i,j)
 839               write (iout,'(2(1pe15.5))')
 840      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 841             enddo
 842             enddo
 843             write (iout,*) 'Lorenz constants'
 844             do k=1,nlor_sccor(i,j)
 845               write (iout,'(3(1pe15.5))')
 846      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 847             enddo
 848           enddo
 849         enddo
 850       endif
 851 C
 852 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 853 C         interaction energy of the Gly, Ala, and Pro prototypes.
 854 C
 855       write (iout,*) "Reading Fourier"
 856       call flush(iout)
 857       read (ifourier,*) nloctyp
 858       do i=0,nloctyp-1
 859         read (ifourier,*)
 860         read (ifourier,*) (b(ii,i),ii=1,13)
 861         if (lprint) then
 862         write (iout,*) 'Type',i
 863         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 864         endif
 865       write (iout,*) "Reading Fourier 1"
 866       call flush(iout)
 867 #ifdef NEWCORR
 868         read (ifourier,*) (bnew1(ii,1,i),ii=1,3)
 869         read (ifourier,*) (bnew2(ii,1,i),ii=1,3)
 870         read (ifourier,*) (bnew1(ii,2,i),ii=1,1)
 871         read (ifourier,*) (bnew2(ii,2,i),ii=1,1)
 872         read (ifourier,*) (eenew(ii,i),ii=1,1)
 873         CC(1,1,i)= b(7,i)
 874         CC(2,2,i)=-b(7,i)
 875         CC(2,1,i)= b(9,i)
 876         CC(1,2,i)= b(9,i)
 877         CC(1,1,-i)= b(7,i)
 878         CC(2,2,-i)=-b(7,i)
 879         CC(2,1,-i)=-b(9,i)
 880         CC(1,2,-i)=-b(9,i)
 881         Ctilde(1,1,i)=b(7,i)
 882         Ctilde(1,2,i)=b(9,i)
 883         Ctilde(2,1,i)=-b(9,i)
 884         Ctilde(2,2,i)=b(7,i)
 885         Ctilde(1,1,-i)=b(7,i)
 886         Ctilde(1,2,-i)=-b(9,i)
 887         Ctilde(2,1,-i)=b(9,i)
 888         Ctilde(2,2,-i)=b(7,i)
 889
 890         DD(1,1,i)= b(6,i)
 891         DD(2,2,i)=-b(6,i)
 892         DD(2,1,i)= b(8,i)
 893         DD(1,2,i)= b(8,i)
 894         DD(1,1,-i)= b(6,i)
 895         DD(2,2,-i)=-b(6,i)
 896         DD(2,1,-i)=-b(8,i)
 897         DD(1,2,-i)=-b(8,i)
 898         Dtilde(1,1,i)=b(6,i)
 899         Dtilde(1,2,i)=b(8,i)
 900         Dtilde(2,1,i)=-b(8,i)
 901         Dtilde(2,2,i)=b(6,i)
 902         Dtilde(1,1,-i)=b(6,i)
 903         Dtilde(1,2,-i)=-b(8,i)
 904         Dtilde(2,1,-i)=b(8,i)
 905         Dtilde(2,2,-i)=b(6,i)
 906
 907         EEold(1,1,i)= b(10,i)+b(11,i)
 908         eeold(1,1,i)= eenew(1,i)
 909         EEold(2,2,i)=-b(10,i)+b(11,i)
 910         EEold(2,1,i)= b(12,i)-b(13,i)
 911         EEold(1,2,i)= b(12,i)+b(13,i)
 912         EEold(1,1,-i)= b(10,i)+b(11,i)
 913         EEold(2,2,-i)=-b(10,i)+b(11,i)
 914         EEold(2,1,-i)=-b(12,i)+b(13,i)
 915         EEold(1,2,-i)=-b(12,i)-b(13,i)
 916 c        write(iout,*) "TU DOCHODZE"
 917       enddo
 918 c      lprint=.true.
 919       if (lprint) then
 920       do i=1,nloctyp
 921         write (iout,*) 'Type',i
 922         write (iout,*) 'BNEW1(1)'
 923         write (iout,*) (BNEW1(ii,1,i),ii=1,3)
 924         write (iout,*) 'BNEW1(2)'
 925         write (iout,*) BNEW1(1,2,i)
 926         write (iout,*) 'BNEW2(1)'
 927         write (iout,*) (BNEW2(ii,1,i),ii=1,3)
 928         write (iout,*) 'BNEW2(2)'
 929         write (iout,*) BNEW2(1,2,i)
 930         write (iout,*) 'CC'
 931         do j=1,2
 932           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 933         enddo
 934         write(iout,*) 'DD'
 935         do j=1,2
 936           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 937         enddo
 938         write(iout,*) 'EE'
 939         do j=1,2
 940           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
 941         enddo
 942       enddo
 943       endif
 944 #else
 945         B1(1,i)  = b(3,i)
 946         B1(2,i)  = b(5,i)
 947         B1(1,-i) = b(3,i)
 948         B1(2,-i) = -b(5,i)
 949 c        b1(1,i)=0.0d0
 950 c        b1(2,i)=0.0d0
 951         B1tilde(1,i) = b(3,i)
 952         B1tilde(2,i) =-b(5,i)
 953         B1tilde(1,-i) =-b(3,i)
 954         B1tilde(2,-i) =b(5,i)
 955 c        b1tilde(1,i)=0.0d0
 956 c        b1tilde(2,i)=0.0d0
 957         B2(1,i)  = b(2,i)
 958         B2(2,i)  = b(4,i)
 959         B2(1,-i)  =b(2,i)
 960         B2(2,-i)  =-b(4,i)
 961
 962 c        b2(1,i)=0.0d0
 963 c        b2(2,i)=0.0d0
 964         CC(1,1,i)= b(7,i)
 965         CC(2,2,i)=-b(7,i)
 966         CC(2,1,i)= b(9,i)
 967         CC(1,2,i)= b(9,i)
 968         CC(1,1,-i)= b(7,i)
 969         CC(2,2,-i)=-b(7,i)
 970         CC(2,1,-i)=-b(9,i)
 971         CC(1,2,-i)=-b(9,i)
 972 c        CC(1,1,i)=0.0d0
 973 c        CC(2,2,i)=0.0d0
 974 c        CC(2,1,i)=0.0d0
 975 c        CC(1,2,i)=0.0d0
 976         Ctilde(1,1,i)=b(7,i)
 977         Ctilde(1,2,i)=b(9,i)
 978         Ctilde(2,1,i)=-b(9,i)
 979         Ctilde(2,2,i)=b(7,i)
 980         Ctilde(1,1,-i)=b(7,i)
 981         Ctilde(1,2,-i)=-b(9,i)
 982         Ctilde(2,1,-i)=b(9,i)
 983         Ctilde(2,2,-i)=b(7,i)
 984
 985 c        Ctilde(1,1,i)=0.0d0
 986 c        Ctilde(1,2,i)=0.0d0
 987 c        Ctilde(2,1,i)=0.0d0
 988 c        Ctilde(2,2,i)=0.0d0
 989         DD(1,1,i)= b(6,i)
 990         DD(2,2,i)=-b(6,i)
 991         DD(2,1,i)= b(8,i)
 992         DD(1,2,i)= b(8,i)
 993         DD(1,1,-i)= b(6,i)
 994         DD(2,2,-i)=-b(6,i)
 995         DD(2,1,-i)=-b(8,i)
 996         DD(1,2,-i)=-b(8,i)
 997 c        DD(1,1,i)=0.0d0
 998 c        DD(2,2,i)=0.0d0
 999 c        DD(2,1,i)=0.0d0
1000 c        DD(1,2,i)=0.0d0
1001         Dtilde(1,1,i)=b(6,i)
1002         Dtilde(1,2,i)=b(8,i)
1003         Dtilde(2,1,i)=-b(8,i)
1004         Dtilde(2,2,i)=b(6,i)
1005         Dtilde(1,1,-i)=b(6,i)
1006         Dtilde(1,2,-i)=-b(8,i)
1007         Dtilde(2,1,-i)=b(8,i)
1008         Dtilde(2,2,-i)=b(6,i)
1009
1010 c        Dtilde(1,1,i)=0.0d0
1011 c        Dtilde(1,2,i)=0.0d0
1012 c        Dtilde(2,1,i)=0.0d0
1013 c        Dtilde(2,2,i)=0.0d0
1014         EE(1,1,i)= b(10,i)+b(11,i)
1015         EE(2,2,i)=-b(10,i)+b(11,i)
1016         EE(2,1,i)= b(12,i)-b(13,i)
1017         EE(1,2,i)= b(12,i)+b(13,i)
1018         EE(1,1,-i)= b(10,i)+b(11,i)
1019         EE(2,2,-i)=-b(10,i)+b(11,i)
1020         EE(2,1,-i)=-b(12,i)+b(13,i)
1021         EE(1,2,-i)=-b(12,i)-b(13,i)
1022
1023 c        ee(1,1,i)=1.0d0
1024 c        ee(2,2,i)=1.0d0
1025 c        ee(2,1,i)=0.0d0
1026 c        ee(1,2,i)=0.0d0
1027 c        ee(2,1,i)=ee(1,2,i)
1028
1029       enddo
1030       if (lprint) then
1031       do i=1,nloctyp
1032         write (iout,*) 'Type',i
1033         write (iout,*) 'B1'
1034 c        write (iout,'(f10.5)') B1(:,i)
1035         write(iout,*) B1(1,i),B1(2,i)
1036         write (iout,*) 'B2'
1037 c        write (iout,'(f10.5)') B2(:,i)
1038         write(iout,*) B2(1,i),B2(2,i)
1039         write (iout,*) 'CC'
1040         do j=1,2
1041           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1042         enddo
1043         write(iout,*) 'DD'
1044         do j=1,2
1045           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1046         enddo
1047         write(iout,*) 'EE'
1048         do j=1,2
1049           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1050         enddo
1051       enddo
1052       endif
1053 #endif
1054 C
1055 C Read electrostatic-interaction parameters
1056 C
1057       if (lprint) then
1058         write (iout,'(/a)') 'Electrostatic interaction constants:'
1059         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1060      &            'IT','JT','APP','BPP','AEL6','AEL3'
1061       endif
1062       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1063       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1064       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1065       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1066       close (ielep)
1067       do i=1,2
1068         do j=1,2
1069         rri=rpp(i,j)**6
1070         app (i,j)=epp(i,j)*rri*rri
1071         bpp (i,j)=-2.0D0*epp(i,j)*rri
1072         ael6(i,j)=elpp6(i,j)*4.2D0**6
1073         ael3(i,j)=elpp3(i,j)*4.2D0**3
1074         lprint=.true.
1075         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1076      &                    ael6(i,j),ael3(i,j)
1077          lprint=.false.
1078         enddo
1079       enddo
1080 C
1081 C Read side-chain interaction parameters.
1082 C
1083       read (isidep,*) ipot,expon
1084       if (ipot.lt.1 .or. ipot.gt.5) then
1085         write (iout,'(2a)') 'Error while reading SC interaction',
1086      &               'potential file - unknown potential type.'
1087         stop
1088       endif
1089       expon2=expon/2
1090       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1091      & ', exponents are ',expon,2*expon
1092       goto (10,20,30,30,40) ipot
1093 C----------------------- LJ potential ---------------------------------
1094    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1095       if (lprint) then
1096         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1097         write (iout,'(a/)') 'The epsilon array:'
1098         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1099         write (iout,'(/a)') 'One-body parameters:'
1100         write (iout,'(a,4x,a)') 'residue','sigma'
1101         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1102       endif
1103       goto 50
1104 C----------------------- LJK potential --------------------------------
1105    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1106      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1107       if (lprint) then
1108         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1109         write (iout,'(a/)') 'The epsilon array:'
1110         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1111         write (iout,'(/a)') 'One-body parameters:'
1112         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1113         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1114      &        i=1,ntyp)
1115       endif
1116       goto 50
1117 C---------------------- GB or BP potential -----------------------------
1118    30 do i=1,ntyp
1119        read (isidep,*)(eps(i,j),j=i,ntyp)
1120       enddo
1121       read (isidep,*)(sigma0(i),i=1,ntyp)
1122       read (isidep,*)(sigii(i),i=1,ntyp)
1123       read (isidep,*)(chip(i),i=1,ntyp)
1124       read (isidep,*)(alp(i),i=1,ntyp)
1125 C For the GB potential convert sigma'**2 into chi'
1126       if (ipot.eq.4) then
1127         do i=1,ntyp
1128           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1129         enddo
1130       endif
1131       if (lprint) then
1132         write (iout,'(/a/)') 'Parameters of the BP potential:'
1133         write (iout,'(a/)') 'The epsilon array:'
1134         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1135         write (iout,'(/a)') 'One-body parameters:'
1136         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1137      &       '    chip  ','    alph  '
1138         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1139      &                     chip(i),alp(i),i=1,ntyp)
1140       endif
1141       goto 50
1142 C--------------------- GBV potential -----------------------------------
1143    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1144      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1145      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1146       if (lprint) then
1147         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1148         write (iout,'(a/)') 'The epsilon array:'
1149         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1150         write (iout,'(/a)') 'One-body parameters:'
1151         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1152      &      's||/s_|_^2','    chip  ','    alph  '
1153         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1154      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1155       endif
1156    50 continue
1157       close (isidep)
1158 C-----------------------------------------------------------------------
1159 C Calculate the "working" parameters of SC interactions.
1160       do i=2,ntyp
1161         do j=1,i-1
1162           eps(i,j)=eps(j,i)
1163         enddo
1164       enddo
1165       do i=1,ntyp
1166         do j=i,ntyp
1167           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1168           sigma(j,i)=sigma(i,j)
1169           rs0(i,j)=dwa16*sigma(i,j)
1170           rs0(j,i)=rs0(i,j)
1171         enddo
1172       enddo
1173       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1174      & 'Working parameters of the SC interactions:',
1175      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1176      & '  chi1   ','   chi2   '
1177       do i=1,ntyp
1178         do j=i,ntyp
1179           epsij=eps(i,j)
1180           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1181             rrij=sigma(i,j)
1182           else
1183             rrij=rr0(i)+rr0(j)
1184           endif
1185           r0(i,j)=rrij
1186           r0(j,i)=rrij
1187           rrij=rrij**expon
1188           epsij=eps(i,j)
1189           sigeps=dsign(1.0D0,epsij)
1190           epsij=dabs(epsij)
1191           aa(i,j)=epsij*rrij*rrij
1192           bb(i,j)=-sigeps*epsij*rrij
1193           aa(j,i)=aa(i,j)
1194           bb(j,i)=bb(i,j)
1195           if (ipot.gt.2) then
1196             sigt1sq=sigma0(i)**2
1197             sigt2sq=sigma0(j)**2
1198             sigii1=sigii(i)
1199             sigii2=sigii(j)
1200             ratsig1=sigt2sq/sigt1sq
1201             ratsig2=1.0D0/ratsig1
1202             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1203             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1204             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1205           else
1206             rsum_max=sigma(i,j)
1207           endif
1208 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1209             sigmaii(i,j)=rsum_max
1210             sigmaii(j,i)=rsum_max
1211 c         else
1212 c           sigmaii(i,j)=r0(i,j)
1213 c           sigmaii(j,i)=r0(i,j)
1214 c         endif
1215 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1216           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1217             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1218             augm(i,j)=epsij*r_augm**(2*expon)
1219 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1220             augm(j,i)=augm(i,j)
1221           else
1222             augm(i,j)=0.0D0
1223             augm(j,i)=0.0D0
1224           endif
1225           if (lprint) then
1226             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1227      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1228      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1229           endif
1230         enddo
1231       enddo
1232 C
1233 C Define the SC-p interaction constants
1234 C
1235 #ifdef OLDSCP
1236       do i=1,20
1237 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1238 C helix formation)
1239 c       aad(i,1)=0.3D0*4.0D0**12
1240 C Following line for constants currently implemented
1241 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1242         aad(i,1)=1.5D0*4.0D0**12
1243 c       aad(i,1)=0.17D0*5.6D0**12
1244         aad(i,2)=aad(i,1)
1245 C "Soft" SC-p repulsion
1246         bad(i,1)=0.0D0
1247 C Following line for constants currently implemented
1248 c       aad(i,1)=0.3D0*4.0D0**6
1249 C "Hard" SC-p repulsion
1250         bad(i,1)=3.0D0*4.0D0**6
1251 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1252         bad(i,2)=bad(i,1)
1253 c       aad(i,1)=0.0D0
1254 c       aad(i,2)=0.0D0
1255 c       bad(i,1)=1228.8D0
1256 c       bad(i,2)=1228.8D0
1257       enddo
1258 #else
1259 C
1260 C 8/9/01 Read the SC-p interaction constants from file
1261 C
1262       do i=1,ntyp
1263         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1264       enddo
1265       do i=1,ntyp
1266         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1267         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1268         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1269         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1270       enddo
1271
1272       if (lprint) then
1273         write (iout,*) "Parameters of SC-p interactions:"
1274         do i=1,20
1275           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1276      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1277         enddo
1278       endif
1279 #endif
1280 C
1281 C Define the constants of the disulfide bridge
1282 C
1283       ebr=-5.50D0
1284 c
1285 c Old arbitrary potential - commented out.
1286 c
1287 c      dbr= 4.20D0
1288 c      fbr= 3.30D0
1289 c
1290 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1291 c energy surface of diethyl disulfide.
1292 c A. Liwo and U. Kozlowska, 11/24/03
1293 c
1294       D0CM = 3.78d0
1295       AKCM = 15.1d0
1296       AKTH = 11.0d0
1297       AKCT = 12.0d0
1298       V1SS =-1.08d0
1299       V2SS = 7.61d0
1300       V3SS = 13.7d0
1301
1302       if (lprint) then
1303       write (iout,'(/a)') "Disulfide bridge parameters:"
1304       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1305       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1306       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1307       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1308      & ' v3ss:',v3ss
1309       endif
1310       return
1311       end