source/wham/src-M/parmread.F

   1       subroutine parmread(iparm,*)
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8       include 'DIMENSIONS.ZSCOPT'
   9       include 'DIMENSIONS.FREE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.CHAIN'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.GEO'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.TORSION'
  16       include 'COMMON.FFIELD'
  17       include 'COMMON.NAMES'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.WEIGHTS'
  20       include 'COMMON.ENEPS'
  21       include 'COMMON.SCCOR'
  22       include 'COMMON.SCROT'
  23       include 'COMMON.FREE'
  24       include 'COMMON.SHIELD'
  25       include 'COMMON.CONTROL'
  26       character*1 t1,t2,t3
  27       character*1 onelett(4) /"G","A","P","D"/
  28       character*1 toronelet(-2:2) /"p","a","G","A","P"/
  29       logical lprint
  30       dimension blower(3,3,maxlob)
  31       character*800 controlcard
  32       character*256 bondname_t,thetname_t,rotname_t,torname_t,
  33      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
  34      & sccorname_t
  35       integer ilen
  36       external ilen
  37       character*16 key
  38       integer iparm
  39       double precision ip,mp
  40       character*6 res1
  41       character*3 lancuch,ucase
  42 C      write (iout,*) "KURWA"
  43 C
  44 C Body
  45 C
  46       call getenv("PRINT_PARM",lancuch)
  47       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  48 C Set LPRINT=.TRUE. for debugging
  49       dwa16=2.0d0**(1.0d0/6.0d0)
  50       itypro=20
  51 C Assign virtual-bond length
  52       vbl=3.8D0
  53       vblinv=1.0D0/vbl
  54       vblinv2=vblinv*vblinv
  55       call card_concat(controlcard,.true.)
  56       wname(4)="WCORRH"
  57       do i=1,n_ene
  58         key = wname(i)(:ilen(wname(i)))
  59         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
  60       enddo
  61
  62       write (iout,*) "iparm",iparm," myparm",myparm
  63 c If reading not own parameters, skip assignment
  64       call reada(controlcard,"D0CM",d0cm,3.78d0)
  65       call reada(controlcard,"AKCM",akcm,15.1d0)
  66       call reada(controlcard,"AKTH",akth,11.0d0)
  67       call reada(controlcard,"AKCT",akct,12.0d0)
  68       call reada(controlcard,"V1SS",v1ss,-1.08d0)
  69       call reada(controlcard,"V2SS",v2ss,7.61d0)
  70       call reada(controlcard,"V3SS",v3ss,13.7d0)
  71       call reada(controlcard,"EBR",ebr,-5.50D0)
  72       call reada(controlcard,"DTRISS",dtriss,1.0D0)
  73       call reada(controlcard,"ATRISS",atriss,0.3D0)
  74       call reada(controlcard,"BTRISS",btriss,0.02D0)
  75       call reada(controlcard,"CTRISS",ctriss,1.0D0)
  76 c      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
  77       write (iout,*) "DYN_SS",dyn_ss
  78       write(iout,*) "ATRISS",atriss
  79       write(iout,*) "BTRISS",btriss
  80       write(iout,*) "CTRISS",ctriss
  81       write(iout,*) "DTRISS",dtriss
  82       call reada(controlcard,"HT",Ht,0.0D0)
  83 c
  84 c Old arbitrary potential - commented out.
  85 c
  86 c      dbr= 4.20D0
  87 c      fbr= 3.30D0
  88 c
  89
  90       if (iparm.eq.myparm .or. .not.separate_parset) then
  91
  92 c
  93 c Setup weights for UNRES
  94 c
  95       wsc=ww(1)
  96       wscp=ww(2)
  97       welec=ww(3)
  98       wcorr=ww(4)
  99       wcorr5=ww(5)
 100       wcorr6=ww(6)
 101       wel_loc=ww(7)
 102       wturn3=ww(8)
 103       wturn4=ww(9)
 104       wturn6=ww(10)
 105       wang=ww(11)
 106       wscloc=ww(12)
 107       wtor=ww(13)
 108       wtor_d=ww(14)
 109       wvdwpp=ww(16)
 110       wbond=ww(18)
 111       wsccor=ww(19)
 112       whpb=ww(15)
 113       wstrain=ww(15)
 114       wliptran=ww(22)
 115       wshield=ww(25)
 116       endif
 117
 118       call card_concat(controlcard,.false.)
 119
 120 c Return if not own parameters
 121
 122       if (iparm.ne.myparm .and. separate_parset) return
 123
 124       call reads(controlcard,"BONDPAR",bondname_t,bondname)
 125       open (ibond,file=bondname_t,status='old')
 126       rewind(ibond)
 127       call reads(controlcard,"THETPAR",thetname_t,thetname)
 128       open (ithep,file=thetname_t,status='old')
 129       rewind(ithep)
 130       call reads(controlcard,"ROTPAR",rotname_t,rotname)
 131       open (irotam,file=rotname_t,status='old')
 132       rewind(irotam)
 133       call reads(controlcard,"TORPAR",torname_t,torname)
 134       open (itorp,file=torname_t,status='old')
 135       rewind(itorp)
 136       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 137       write (iout,*) "tor_mode",tor_mode
 138       call flush(iout)
 139       if (tor_mode.eq.0)
 140      & open (itordp,file=tordname_t,status='old')
 141       rewind(itordp)
 142       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
 143       open (isccor,file=sccorname_t,status='old')
 144       rewind(isccor)
 145       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 146       open (ifourier,file=fouriername_t,status='old')
 147       rewind(ifourier)
 148       call reads(controlcard,"ELEPAR",elename_t,elename)
 149       open (ielep,file=elename_t,status='old')
 150       rewind(ielep)
 151       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 152       open (isidep,file=sidename_t,status='old')
 153       rewind(isidep)
 154       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 155       open (iscpp,file=scpname_t,status='old')
 156       rewind(iscpp)
 157       write (iout,*) "Parameter set:",iparm
 158       write (iout,*) "Energy-term weights:"
 159       do i=1,n_ene
 160         write (iout,'(a16,f10.5)') wname(i),ww(i)
 161       enddo
 162       write (iout,*) "Sidechain potential file        : ",
 163      &  sidename_t(:ilen(sidename_t))
 164 #ifndef OLDSCP
 165       write (iout,*) "SCp potential file              : ",
 166      &  scpname_t(:ilen(scpname_t))
 167 #endif
 168       write (iout,*) "Electrostatic potential file    : ",
 169      &  elename_t(:ilen(elename_t))
 170       write (iout,*) "Cumulant coefficient file       : ",
 171      &  fouriername_t(:ilen(fouriername_t))
 172       write (iout,*) "Torsional parameter file        : ",
 173      &  torname_t(:ilen(torname_t))
 174       write (iout,*) "Double torsional parameter file : ",
 175      &  tordname_t(:ilen(tordname_t))
 176       write (iout,*) "Backbone-rotamer parameter file : ",
 177      &  sccorname(:ilen(sccorname))
 178       write (iout,*) "Bond & inertia constant file    : ",
 179      &  bondname_t(:ilen(bondname_t))
 180       write (iout,*) "Bending parameter file          : ",
 181      &  thetname_t(:ilen(thetname_t))
 182       write (iout,*) "Rotamer parameter file          : ",
 183      &  rotname_t(:ilen(rotname_t))
 184
 185 c
 186 c Read the virtual-bond parameters, masses, and moments of inertia
 187 c and Stokes' radii of the peptide group and side chains
 188 c
 189 #ifdef CRYST_BOND
 190       read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp
 191       do i=1,ntyp
 192         nbondterm(i)=1
 193         read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i)
 194         dsc(i) = vbldsc0(1,i)
 195         if (i.eq.10) then
 196           dsc_inv(i)=0.0D0
 197         else
 198           dsc_inv(i)=1.0D0/dsc(i)
 199         endif
 200       enddo
 201 #else
 202       read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk
 203       do i=1,ntyp
 204         read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
 205      &   aksc(j,i),abond0(j,i),j=1,nbondterm(i))
 206         dsc(i) = vbldsc0(1,i)
 207         if (i.eq.10) then
 208           dsc_inv(i)=0.0D0
 209         else
 210           dsc_inv(i)=1.0D0/dsc(i)
 211         endif
 212       enddo
 213 #endif
 214       if (lprint) then
 215         write(iout,'(/a/)')"Force constants virtual bonds:"
 216         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
 217      &   'inertia','Pstok'
 218         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 219         do i=1,ntyp
 220           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
 221      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
 222           do j=2,nbondterm(i)
 223             write (iout,'(13x,3f10.5)')
 224      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
 225           enddo
 226         enddo
 227       endif
 228        write (iout,*) "iliptranpar",iliptranpar
 229        write (iout,*) "liptranname ",liptranname
 230        read(iliptranpar,*,end=1161,err=1161) pepliptran
 231        write (iout,*) "pepliptran",pepliptran
 232        do i=1,ntyp
 233        read(iliptranpar,*,end=1161,err=1161) liptranene(i)
 234        write (iout,*) i,liptranene(i)
 235        enddo
 236        rewind iliptranpar
 237 #ifdef CRYST_THETA
 238 C
 239 C Read the parameters of the probability distribution/energy expression
 240 C of the virtual-bond valence angles theta
 241 C
 242       do i=1,ntyp
 243         read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
 244      &    (bthet(j,i,1,1),j=1,2)
 245         read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
 246         read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
 247         read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
 248         sigc0(i)=sigc0(i)**2
 249       enddo
 250       do i=1,ntyp
 251       athet(1,i,1,-1)=athet(1,i,1,1)
 252       athet(2,i,1,-1)=athet(2,i,1,1)
 253       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 254       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 255       athet(1,i,-1,1)=-athet(1,i,1,1)
 256       athet(2,i,-1,1)=-athet(2,i,1,1)
 257       bthet(1,i,-1,1)=bthet(1,i,1,1)
 258       bthet(2,i,-1,1)=bthet(2,i,1,1)
 259       enddo
 260       do i=-ntyp,-1
 261       a0thet(i)=a0thet(-i)
 262       athet(1,i,-1,-1)=athet(1,-i,1,1)
 263       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 264       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 265       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 266       athet(1,i,-1,1)=athet(1,-i,1,1)
 267       athet(2,i,-1,1)=-athet(2,-i,1,1)
 268       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 269       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 270       athet(1,i,1,-1)=-athet(1,-i,1,1)
 271       athet(2,i,1,-1)=athet(2,-i,1,1)
 272       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 273       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 274       theta0(i)=theta0(-i)
 275       sig0(i)=sig0(-i)
 276       sigc0(i)=sigc0(-i)
 277        do j=0,3
 278         polthet(j,i)=polthet(j,-i)
 279        enddo
 280        do j=1,3
 281          gthet(j,i)=gthet(j,-i)
 282        enddo
 283       enddo
 284       close (ithep)
 285       if (lprint) then
 286 c       write (iout,'(a)')
 287 c    &    'Parameters of the virtual-bond valence angles:'
 288 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 289 c    & '    ATHETA0   ','         A1   ','        A2    ',
 290 c    & '        B1    ','         B2   '
 291 c       do i=1,ntyp
 292 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 293 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 294 c       enddo
 295 c       write (iout,'(/a/9x,5a/79(1h-))')
 296 c    & 'Parameters of the expression for sigma(theta_c):',
 297 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 298 c    & '     ALPH3    ','    SIGMA0C   '
 299 c       do i=1,ntyp
 300 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 301 c    &      (polthet(j,i),j=0,3),sigc0(i)
 302 c       enddo
 303 c       write (iout,'(/a/9x,5a/79(1h-))')
 304 c    & 'Parameters of the second gaussian:',
 305 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 306 c    & '        G2    ','         G3   '
 307 c       do i=1,ntyp
 308 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 309 c    &       sig0(i),(gthet(j,i),j=1,3)
 310 c       enddo
 311         write (iout,'(a)')
 312      &    'Parameters of the virtual-bond valence angles:'
 313         write (iout,'(/a/9x,5a/79(1h-))')
 314      & 'Coefficients of expansion',
 315      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 316      & '   b1*10^1    ','    b2*10^1   '
 317         do i=1,ntyp
 318           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 319      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
 320      &        (10*bthet(j,i,1,1),j=1,2)
 321         enddo
 322         write (iout,'(/a/9x,5a/79(1h-))')
 323      & 'Parameters of the expression for sigma(theta_c):',
 324      & ' alpha0       ','  alph1       ',' alph2        ',
 325      & ' alhp3        ','   sigma0c    '
 326         do i=1,ntyp
 327           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 328      &      (polthet(j,i),j=0,3),sigc0(i)
 329         enddo
 330         write (iout,'(/a/9x,5a/79(1h-))')
 331      & 'Parameters of the second gaussian:',
 332      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 333      & '        G2    ','   G3*10^1    '
 334         do i=1,ntyp
 335           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 336      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 337         enddo
 338       endif
 339 #else
 340       IF (tor_mode.eq.0) THEN
 341 C
 342 C Read the parameters of Utheta determined from ab initio surfaces
 343 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 344 C
 345       read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
 346      &  ntheterm3,nsingle,ndouble
 347       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 348       read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
 349       do i=-ntyp1,-1
 350         ithetyp(i)=-ithetyp(-i)
 351       enddo
 352       do iblock=1,2
 353       do i=-maxthetyp,maxthetyp
 354         do j=-maxthetyp,maxthetyp
 355           do k=-maxthetyp,maxthetyp
 356             aa0thet(i,j,k,iblock)=0.0d0
 357             do l=1,ntheterm
 358               aathet(l,i,j,k,iblock)=0.0d0
 359             enddo
 360             do l=1,ntheterm2
 361               do m=1,nsingle
 362                 bbthet(m,l,i,j,k,iblock)=0.0d0
 363                 ccthet(m,l,i,j,k,iblock)=0.0d0
 364                 ddthet(m,l,i,j,k,iblock)=0.0d0
 365                 eethet(m,l,i,j,k,iblock)=0.0d0
 366               enddo
 367             enddo
 368             do l=1,ntheterm3
 369               do m=1,ndouble
 370                 do mm=1,ndouble
 371                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 372                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 373                 enddo
 374               enddo
 375             enddo
 376           enddo
 377         enddo
 378       enddo
 379       enddo
 380 C      write (iout,*) "KURWA1"
 381       do iblock=1,2
 382       do i=0,nthetyp
 383         do j=-nthetyp,nthetyp
 384           do k=-nthetyp,nthetyp
 385             read (ithep,'(6a)',end=111,err=111) res1
 386             write(iout,*) res1,i,j,k
 387             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
 388             read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
 389      &        l=1,ntheterm)
 390             read (ithep,*,end=111,err=111)
 391      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 392      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 393      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 394      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 395      &        ,ll=1,ntheterm2)
 396             read (ithep,*,end=111,err=111)
 397      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 398      &      ffthet(lll,llll,ll,i,j,k,iblock),
 399      &         ggthet(llll,lll,ll,i,j,k,iblock)
 400      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 401      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 402           enddo
 403         enddo
 404       enddo
 405 C       write(iout,*) "KURWA1.1"
 406 C
 407 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 408 C coefficients of theta-and-gamma-dependent terms are zero.
 409 C
 410       do i=1,nthetyp
 411         do j=1,nthetyp
 412           do l=1,ntheterm
 413             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 414             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 415           enddo
 416           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 417           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 418         enddo
 419         do l=1,ntheterm
 420           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 421         enddo
 422         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 423       enddo
 424       enddo
 425 C       write(iout,*) "KURWA1.5"
 426 C Substitution for D aminoacids from symmetry.
 427       do iblock=1,2
 428       do i=-nthetyp,0
 429         do j=-nthetyp,nthetyp
 430           do k=-nthetyp,nthetyp
 431            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 432            do l=1,ntheterm
 433            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 434            enddo
 435            do ll=1,ntheterm2
 436             do lll=1,nsingle
 437             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 438             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 439             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 440             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 441             enddo
 442           enddo
 443           do ll=1,ntheterm3
 444            do lll=2,ndouble
 445             do llll=1,lll-1
 446             ffthet(llll,lll,ll,i,j,k,iblock)=
 447      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 448             ffthet(lll,llll,ll,i,j,k,iblock)=
 449      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 450             ggthet(llll,lll,ll,i,j,k,iblock)=
 451      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 452             ggthet(lll,llll,ll,i,j,k,iblock)=
 453      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 454             enddo !ll
 455            enddo  !lll
 456           enddo   !llll
 457          enddo    !k
 458         enddo     !j
 459        enddo      !i
 460       enddo       !iblock
 461
 462 C
 463 C Control printout of the coefficients of virtual-bond-angle potentials
 464 C
 465       if (lprint) then
 466         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 467         do iblock=1,2
 468         do i=1,nthetyp+1
 469           do j=1,nthetyp+1
 470             do k=1,nthetyp+1
 471               write (iout,'(//4a)')
 472      &         'Type ',onelett(i),onelett(j),onelett(k)
 473               write (iout,'(//a,10x,a)') " l","a[l]"
 474               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 475               write (iout,'(i2,1pe15.5)')
 476      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 477             do l=1,ntheterm2
 478               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 479      &          "b",l,"c",l,"d",l,"e",l
 480               do m=1,nsingle
 481                 write (iout,'(i2,4(1pe15.5))') m,
 482      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 483      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 484               enddo
 485             enddo
 486             do l=1,ntheterm3
 487               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 488      &          "f+",l,"f-",l,"g+",l,"g-",l
 489               do m=2,ndouble
 490                 do n=1,m-1
 491                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 492      &              ffthet(n,m,l,i,j,k,iblock),
 493      &              ffthet(m,n,l,i,j,k,iblock),
 494      &              ggthet(n,m,l,i,j,k,iblock),
 495      &              ggthet(m,n,l,i,j,k,iblock)
 496                 enddo
 497               enddo
 498             enddo
 499           enddo
 500         enddo
 501       enddo
 502       enddo
 503       call flush(iout)
 504       endif
 505
 506       ELSE
 507
 508 C here will be the apropriate recalibrating for D-aminoacid
 509       read (ithep,*,end=111,err=111) nthetyp
 510       do i=-nthetyp+1,nthetyp-1
 511         read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
 512         do j=0,nbend_kcc_Tb(i)
 513           read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
 514         enddo
 515       enddo
 516       if (lprint) then
 517         write (iout,'(a)')
 518      &    "Parameters of the valence-only potentials"
 519         do i=-nthetyp+1,nthetyp-1
 520         write (iout,'(2a)') "Type ",toronelet(i)
 521         do k=0,nbend_kcc_Tb(i)
 522           write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
 523         enddo
 524         enddo
 525       endif
 526
 527       ENDIF ! TOR_MODE
 528
 529       close(ithep)
 530 #endif
 531 C      write(iout,*) 'KURWA2'
 532 #ifdef CRYST_SC
 533 C
 534 C Read the parameters of the probability distribution/energy expression
 535 C of the side chains.
 536 C
 537       do i=1,ntyp
 538 cc      write (iout,*) "tu dochodze",i
 539         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 540         if (i.eq.10) then
 541           dsc_inv(i)=0.0D0
 542         else
 543           dsc_inv(i)=1.0D0/dsc(i)
 544         endif
 545         if (i.ne.10) then
 546         do j=1,nlob(i)
 547           do k=1,3
 548             do l=1,3
 549               blower(l,k,j)=0.0D0
 550             enddo
 551           enddo
 552         enddo
 553         bsc(1,i)=0.0D0
 554         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 555      &    ((blower(k,l,1),l=1,k),k=1,3)
 556         censc(1,1,-i)=censc(1,1,i)
 557         censc(2,1,-i)=censc(2,1,i)
 558         censc(3,1,-i)=-censc(3,1,i)
 559         do j=2,nlob(i)
 560           read (irotam,*,end=112,err=112) bsc(j,i)
 561           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 562      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 563         censc(1,j,-i)=censc(1,j,i)
 564         censc(2,j,-i)=censc(2,j,i)
 565         censc(3,j,-i)=-censc(3,j,i)
 566 C BSC is amplitude of Gaussian
 567         enddo
 568         do j=1,nlob(i)
 569           do k=1,3
 570             do l=1,k
 571               akl=0.0D0
 572               do m=1,3
 573                 akl=akl+blower(k,m,j)*blower(l,m,j)
 574               enddo
 575               gaussc(k,l,j,i)=akl
 576               gaussc(l,k,j,i)=akl
 577              if (((k.eq.3).and.(l.ne.3))
 578      &        .or.((l.eq.3).and.(k.ne.3))) then
 579                 gaussc(k,l,j,-i)=-akl
 580                 gaussc(l,k,j,-i)=-akl
 581               else
 582                 gaussc(k,l,j,-i)=akl
 583                 gaussc(l,k,j,-i)=akl
 584               endif
 585             enddo
 586           enddo
 587         enddo
 588         endif
 589       enddo
 590       close (irotam)
 591       if (lprint) then
 592         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 593         do i=1,ntyp
 594           nlobi=nlob(i)
 595           if (nlobi.gt.0) then
 596           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 597      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 598 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 599 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 600 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 601 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 602            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 603      &                             'log h',(bsc(j,i),j=1,nlobi)
 604            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 605      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 606 c          write (iout,'(a)')
 607 c         do j=1,nlobi
 608 c           ind=0
 609 c           do k=1,3
 610 c             do l=1,k
 611 c              ind=ind+1
 612 c              blower(k,l,j)=gaussc(ind,j,i)
 613 c             enddo
 614 c           enddo
 615 c         enddo
 616           do k=1,3
 617             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 618      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 619           enddo
 620           endif
 621         enddo
 622       endif
 623 #else
 624 C
 625 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 626 C added by Urszula Kozlowska 07/11/2007
 627 C
 628       do i=1,ntyp
 629         read (irotam,*,end=112,err=112)
 630        if (i.eq.10) then
 631          read (irotam,*,end=112,err=112)
 632        else
 633          do j=1,65
 634            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 635          enddo
 636        endif
 637       enddo
 638 #endif
 639       close(irotam)
 640 C
 641 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 642 C         interaction energy of the Gly, Ala, and Pro prototypes.
 643 C
 644       read (ifourier,*,end=115,err=115) nloctyp
 645       SPLIT_FOURIERTOR = nloctyp.lt.0
 646       nloctyp = iabs(nloctyp)
 647 #ifdef NEWCORR
 648       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
 649       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
 650       itype2loc(ntyp1)=nloctyp
 651       iloctyp(nloctyp)=ntyp1
 652       do i=1,ntyp1
 653         itype2loc(-i)=-itype2loc(i)
 654       enddo
 655 #else
 656       iloctyp(0)=10
 657       iloctyp(1)=9
 658       iloctyp(2)=20
 659       iloctyp(3)=ntyp1
 660 #endif
 661       do i=1,nloctyp
 662         iloctyp(-i)=-iloctyp(i)
 663       enddo
 664 c      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
 665 c      write (iout,*) "nloctyp",nloctyp,
 666 c     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
 667 #ifdef NEWCORR
 668       do i=0,nloctyp-1
 669 c             write (iout,*) "NEWCORR",i
 670         read (ifourier,*,end=115,err=115)
 671         do ii=1,3
 672           do j=1,2
 673             read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
 674           enddo
 675         enddo
 676 c             write (iout,*) "NEWCORR BNEW1"
 677 c             write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
 678         do ii=1,3
 679           do j=1,2
 680             read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
 681           enddo
 682         enddo
 683 c             write (iout,*) "NEWCORR BNEW2"
 684 c             write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
 685         do kk=1,3
 686           read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
 687           read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
 688         enddo
 689 c             write (iout,*) "NEWCORR CCNEW"
 690 c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
 691         do kk=1,3
 692           read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
 693           read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
 694         enddo
 695 c             write (iout,*) "NEWCORR DDNEW"
 696 c             write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
 697         do ii=1,2
 698           do jj=1,2
 699             do kk=1,2
 700               read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
 701             enddo
 702           enddo
 703         enddo
 704 c             write (iout,*) "NEWCORR EENEW1"
 705 c             write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
 706         do ii=1,3
 707           read (ifourier,*,end=115,err=115) e0new(ii,i)
 708         enddo
 709 c             write (iout,*) (e0new(ii,i),ii=1,3)
 710       enddo
 711 c             write (iout,*) "NEWCORR EENEW"
 712       do i=0,nloctyp-1
 713         do ii=1,3
 714           ccnew(ii,1,i)=ccnew(ii,1,i)/2
 715           ccnew(ii,2,i)=ccnew(ii,2,i)/2
 716           ddnew(ii,1,i)=ddnew(ii,1,i)/2
 717           ddnew(ii,2,i)=ddnew(ii,2,i)/2
 718         enddo
 719       enddo
 720       do i=1,nloctyp-1
 721         do ii=1,3
 722           bnew1(ii,1,-i)= bnew1(ii,1,i)
 723           bnew1(ii,2,-i)=-bnew1(ii,2,i)
 724           bnew2(ii,1,-i)= bnew2(ii,1,i)
 725           bnew2(ii,2,-i)=-bnew2(ii,2,i)
 726         enddo
 727         do ii=1,3
 728 c          ccnew(ii,1,i)=ccnew(ii,1,i)/2
 729 c          ccnew(ii,2,i)=ccnew(ii,2,i)/2
 730 c          ddnew(ii,1,i)=ddnew(ii,1,i)/2
 731 c          ddnew(ii,2,i)=ddnew(ii,2,i)/2
 732           ccnew(ii,1,-i)=ccnew(ii,1,i)
 733           ccnew(ii,2,-i)=-ccnew(ii,2,i)
 734           ddnew(ii,1,-i)=ddnew(ii,1,i)
 735           ddnew(ii,2,-i)=-ddnew(ii,2,i)
 736         enddo
 737         e0new(1,-i)= e0new(1,i)
 738         e0new(2,-i)=-e0new(2,i)
 739         e0new(3,-i)=-e0new(3,i)
 740         do kk=1,2
 741           eenew(kk,1,1,-i)= eenew(kk,1,1,i)
 742           eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
 743           eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
 744           eenew(kk,2,2,-i)= eenew(kk,2,2,i)
 745         enddo
 746       enddo
 747       if (lprint) then
 748         write (iout,'(a)') "Coefficients of the multibody terms"
 749         do i=-nloctyp+1,nloctyp-1
 750           write (iout,*) "Type: ",onelet(iloctyp(i))
 751           write (iout,*) "Coefficients of the expansion of B1"
 752           do j=1,2
 753             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
 754           enddo
 755           write (iout,*) "Coefficients of the expansion of B2"
 756           do j=1,2
 757             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
 758           enddo
 759           write (iout,*) "Coefficients of the expansion of C"
 760           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
 761           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
 762           write (iout,*) "Coefficients of the expansion of D"
 763           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
 764           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
 765           write (iout,*) "Coefficients of the expansion of E"
 766           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
 767           do j=1,2
 768             do k=1,2
 769               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
 770             enddo
 771           enddo
 772         enddo
 773       endif
 774       IF (SPLIT_FOURIERTOR) THEN
 775       do i=0,nloctyp-1
 776 c             write (iout,*) "NEWCORR TOR",i
 777         read (ifourier,*,end=115,err=115)
 778         do ii=1,3
 779           do j=1,2
 780             read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
 781           enddo
 782         enddo
 783 c             write (iout,*) "NEWCORR BNEW1 TOR"
 784 c             write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
 785         do ii=1,3
 786           do j=1,2
 787             read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
 788           enddo
 789         enddo
 790 c             write (iout,*) "NEWCORR BNEW2 TOR"
 791 c             write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
 792         do kk=1,3
 793           read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
 794           read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
 795         enddo
 796 c             write (iout,*) "NEWCORR CCNEW TOR"
 797 c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
 798         do kk=1,3
 799           read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
 800           read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
 801         enddo
 802 c             write (iout,*) "NEWCORR DDNEW TOR"
 803 c             write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
 804         do ii=1,2
 805           do jj=1,2
 806             do kk=1,2
 807               read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
 808             enddo
 809           enddo
 810         enddo
 811 c         write (iout,*) "NEWCORR EENEW1 TOR"
 812 c         write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
 813         do ii=1,3
 814           read (ifourier,*,end=115,err=115) e0newtor(ii,i)
 815         enddo
 816 c             write (iout,*) (e0newtor(ii,i),ii=1,3)
 817       enddo
 818 c             write (iout,*) "NEWCORR EENEW TOR"
 819       do i=0,nloctyp-1
 820         do ii=1,3
 821           ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
 822           ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
 823           ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
 824           ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
 825         enddo
 826       enddo
 827       do i=1,nloctyp-1
 828         do ii=1,3
 829           bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
 830           bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
 831           bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
 832           bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
 833         enddo
 834         do ii=1,3
 835           ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
 836           ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
 837           ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
 838           ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
 839         enddo
 840         e0newtor(1,-i)= e0newtor(1,i)
 841         e0newtor(2,-i)=-e0newtor(2,i)
 842         e0newtor(3,-i)=-e0newtor(3,i)
 843         do kk=1,2
 844           eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
 845           eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
 846           eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
 847           eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
 848         enddo
 849       enddo
 850       if (lprint) then
 851         write (iout,'(a)')
 852      &    "Single-body coefficients of the torsional potentials"
 853         do i=-nloctyp+1,nloctyp-1
 854           write (iout,*) "Type: ",onelet(iloctyp(i))
 855           write (iout,*) "Coefficients of the expansion of B1tor"
 856           do j=1,2
 857             write (iout,'(3hB1(,i1,1h),3f10.5)')
 858      &        j,(bnew1tor(k,j,i),k=1,3)
 859           enddo
 860           write (iout,*) "Coefficients of the expansion of B2tor"
 861           do j=1,2
 862             write (iout,'(3hB2(,i1,1h),3f10.5)')
 863      &        j,(bnew2tor(k,j,i),k=1,3)
 864           enddo
 865           write (iout,*) "Coefficients of the expansion of Ctor"
 866           write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
 867           write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
 868           write (iout,*) "Coefficients of the expansion of Dtor"
 869           write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
 870           write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
 871           write (iout,*) "Coefficients of the expansion of Etor"
 872           write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
 873           do j=1,2
 874             do k=1,2
 875               write (iout,'(1hE,2i1,2f10.5)')
 876      &          j,k,(eenewtor(l,j,k,i),l=1,2)
 877             enddo
 878           enddo
 879         enddo
 880       endif
 881       ELSE
 882       do i=-nloctyp+1,nloctyp-1
 883         do ii=1,3
 884           do j=1,2
 885             bnew1tor(ii,j,i)=bnew1(ii,j,i)
 886           enddo
 887         enddo
 888         do ii=1,3
 889           do j=1,2
 890             bnew2tor(ii,j,i)=bnew2(ii,j,i)
 891           enddo
 892         enddo
 893         do ii=1,3
 894           ccnewtor(ii,1,i)=ccnew(ii,1,i)
 895           ccnewtor(ii,2,i)=ccnew(ii,2,i)
 896           ddnewtor(ii,1,i)=ddnew(ii,1,i)
 897           ddnewtor(ii,2,i)=ddnew(ii,2,i)
 898         enddo
 899       enddo
 900       ENDIF !SPLIT_FOURIER_TOR
 901 #else
 902       if (lprint)
 903      &  write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
 904       do i=0,nloctyp-1
 905         read (ifourier,*,end=115,err=115)
 906         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
 907         if (lprint) then
 908         write (iout,*) 'Type ',onelet(iloctyp(i))
 909         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 910         endif
 911         if (i.gt.0) then
 912         b(2,-i)= b(2,i)
 913         b(3,-i)= b(3,i)
 914         b(4,-i)=-b(4,i)
 915         b(5,-i)=-b(5,i)
 916         endif
 917 c        B1(1,i)  = b(3)
 918 c        B1(2,i)  = b(5)
 919 c        B1(1,-i) = b(3)
 920 c        B1(2,-i) = -b(5)
 921 c        b1(1,i)=0.0d0
 922 c        b1(2,i)=0.0d0
 923 c        B1tilde(1,i) = b(3)
 924 c        B1tilde(2,i) =-b(5)
 925 c        B1tilde(1,-i) =-b(3)
 926 c        B1tilde(2,-i) =b(5)
 927 c        b1tilde(1,i)=0.0d0
 928 c        b1tilde(2,i)=0.0d0
 929 c        B2(1,i)  = b(2)
 930 c        B2(2,i)  = b(4)
 931 c        B2(1,-i)  =b(2)
 932 c        B2(2,-i)  =-b(4)
 933 cc        B1tilde(1,i) = b(3,i)
 934 cc        B1tilde(2,i) =-b(5,i)
 935 C        B1tilde(1,-i) =-b(3,i)
 936 C        B1tilde(2,-i) =b(5,i)
 937 cc        b1tilde(1,i)=0.0d0
 938 cc        b1tilde(2,i)=0.0d0
 939 cc        B2(1,i)  = b(2,i)
 940 cc        B2(2,i)  = b(4,i)
 941 C        B2(1,-i)  =b(2,i)
 942 C        B2(2,-i)  =-b(4,i)
 943
 944 c        b2(1,i)=0.0d0
 945 c        b2(2,i)=0.0d0
 946         CCold(1,1,i)= b(7,i)
 947         CCold(2,2,i)=-b(7,i)
 948         CCold(2,1,i)= b(9,i)
 949         CCold(1,2,i)= b(9,i)
 950         CCold(1,1,-i)= b(7,i)
 951         CCold(2,2,-i)=-b(7,i)
 952         CCold(2,1,-i)=-b(9,i)
 953         CCold(1,2,-i)=-b(9,i)
 954 c        CC(1,1,i)=0.0d0
 955 c        CC(2,2,i)=0.0d0
 956 c        CC(2,1,i)=0.0d0
 957 c        CC(1,2,i)=0.0d0
 958 c        Ctilde(1,1,i)= CCold(1,1,i)
 959 c        Ctilde(1,2,i)= CCold(1,2,i)
 960 c        Ctilde(2,1,i)=-CCold(2,1,i)
 961 c        Ctilde(2,2,i)=-CCold(2,2,i)
 962
 963 c        Ctilde(1,1,i)=0.0d0
 964 c        Ctilde(1,2,i)=0.0d0
 965 c        Ctilde(2,1,i)=0.0d0
 966 c        Ctilde(2,2,i)=0.0d0
 967         DDold(1,1,i)= b(6,i)
 968         DDold(2,2,i)=-b(6,i)
 969         DDold(2,1,i)= b(8,i)
 970         DDold(1,2,i)= b(8,i)
 971         DDold(1,1,-i)= b(6,i)
 972         DDold(2,2,-i)=-b(6,i)
 973         DDold(2,1,-i)=-b(8,i)
 974         DDold(1,2,-i)=-b(8,i)
 975 c        DD(1,1,i)=0.0d0
 976 c        DD(2,2,i)=0.0d0
 977 c        DD(2,1,i)=0.0d0
 978 c        DD(1,2,i)=0.0d0
 979 c        Dtilde(1,1,i)= DD(1,1,i)
 980 c        Dtilde(1,2,i)= DD(1,2,i)
 981 c        Dtilde(2,1,i)=-DD(2,1,i)
 982 c        Dtilde(2,2,i)=-DD(2,2,i)
 983
 984 c        Dtilde(1,1,i)=0.0d0
 985 c        Dtilde(1,2,i)=0.0d0
 986 c        Dtilde(2,1,i)=0.0d0
 987 c        Dtilde(2,2,i)=0.0d0
 988         EEold(1,1,i)= b(10,i)+b(11,i)
 989         EEold(2,2,i)=-b(10,i)+b(11,i)
 990         EEold(2,1,i)= b(12,i)-b(13,i)
 991         EEold(1,2,i)= b(12,i)+b(13,i)
 992         EEold(1,1,-i)= b(10,i)+b(11,i)
 993         EEold(2,2,-i)=-b(10,i)+b(11,i)
 994         EEold(2,1,-i)=-b(12,i)+b(13,i)
 995         EEold(1,2,-i)=-b(12,i)-b(13,i)
 996         write(iout,*) "TU DOCHODZE"
 997         print *,"JESTEM"
 998 c        ee(1,1,i)=1.0d0
 999 c        ee(2,2,i)=1.0d0
1000 c        ee(2,1,i)=0.0d0
1001 c        ee(1,2,i)=0.0d0
1002 c        ee(2,1,i)=ee(1,2,i)
1003       enddo
1004       if (lprint) then
1005       write (iout,*)
1006       write (iout,*)
1007      &"Coefficients of the cumulants (independent of valence angles)"
1008       do i=-nloctyp+1,nloctyp-1
1009         write (iout,*) 'Type ',onelet(iloctyp(i))
1010         write (iout,*) 'B1'
1011         write(iout,'(2f10.5)') B(3,i),B(5,i)
1012         write (iout,*) 'B2'
1013         write(iout,'(2f10.5)') B(2,i),B(4,i)
1014         write (iout,*) 'CC'
1015         do j=1,2
1016           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1017         enddo
1018         write(iout,*) 'DD'
1019         do j=1,2
1020           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1021         enddo
1022         write(iout,*) 'EE'
1023         do j=1,2
1024           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1025         enddo
1026       enddo
1027       endif
1028 #endif
1029 C      write (iout,*) 'KURWAKURWA'
1030 #ifdef CRYST_TOR
1031 C
1032 C Read torsional parameters in old format
1033 C
1034       read (itorp,*,end=113,err=113) ntortyp,nterm_old
1035       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1036       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1037       do i=1,ntortyp
1038         do j=1,ntortyp
1039           read (itorp,'(a)',end=113,err=113)
1040           do k=1,nterm_old
1041             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1042           enddo
1043         enddo
1044       enddo
1045       close (itorp)
1046       if (lprint) then
1047         write (iout,'(/a/)') 'Torsional constants:'
1048         do i=1,ntortyp
1049           do j=1,ntortyp
1050             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1051             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1052           enddo
1053         enddo
1054       endif
1055 #else
1056 C
1057 C Read torsional parameters
1058 C
1059       IF (TOR_MODE.eq.0) THEN
1060
1061       read (itorp,*,end=113,err=113) ntortyp
1062       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1063       do i=1,ntyp1
1064         itype2loc(i)=itortyp(i)
1065       enddo
1066       write (iout,*) 'ntortyp',ntortyp
1067       do i=1,ntyp1
1068         itype2loc(-i)=-itype2loc(i)
1069       enddo
1070       itortyp(ntyp1)=ntortyp
1071       do iblock=1,2
1072       do i=-ntyp,-1
1073        itortyp(i)=-itortyp(-i)
1074       enddo
1075 c      write (iout,*) 'ntortyp',ntortyp
1076       do i=0,ntortyp-1
1077         do j=-ntortyp+1,ntortyp-1
1078           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1079      &          nlor(i,j,iblock)
1080           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1081           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1082           v0ij=0.0d0
1083           si=-1.0d0
1084           do k=1,nterm(i,j,iblock)
1085             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1086      &      v2(k,i,j,iblock)
1087             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1088             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1089             v0ij=v0ij+si*v1(k,i,j,iblock)
1090             si=-si
1091          enddo
1092           do k=1,nlor(i,j,iblock)
1093             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1094      &        vlor2(k,i,j),vlor3(k,i,j)
1095             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1096           enddo
1097           v0(i,j,iblock)=v0ij
1098           v0(-i,-j,iblock)=v0ij
1099         enddo
1100       enddo
1101       enddo
1102       close (itorp)
1103       if (lprint) then
1104         write (iout,'(/a/)') 'Torsional constants:'
1105         do i=1,ntortyp
1106           do j=1,ntortyp
1107             write (iout,*) 'ityp',i,' jtyp',j
1108             write (iout,*) 'Fourier constants'
1109             do k=1,nterm(i,j,iblock)
1110               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1111      &        v2(k,i,j,iblock)
1112             enddo
1113             write (iout,*) 'Lorenz constants'
1114             do k=1,nlor(i,j,iblock)
1115               write (iout,'(3(1pe15.5))')
1116      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1117             enddo
1118           enddo
1119         enddo
1120       endif
1121 C
1122 C 6/23/01 Read parameters for double torsionals
1123 C
1124       do iblock=1,2
1125       do i=0,ntortyp-1
1126         do j=-ntortyp+1,ntortyp-1
1127           do k=-ntortyp+1,ntortyp-1
1128             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1129 c              write (iout,*) "OK onelett",
1130 c     &         i,j,k,t1,t2,t3
1131
1132             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1133      &        .or. t3.ne.toronelet(k)) then
1134               write (iout,*) "Error in double torsional parameter file",
1135      &         i,j,k,t1,t2,t3
1136 #ifdef MPI
1137               call MPI_Finalize(Ierror)
1138 #endif
1139                stop "Error in double torsional parameter file"
1140             endif
1141           read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1142      &         ntermd_2(i,j,k,iblock)
1143             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1144             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1145             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1146      &         ntermd_1(i,j,k,iblock))
1147             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1148      &         ntermd_1(i,j,k,iblock))
1149             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1150      &         ntermd_1(i,j,k,iblock))
1151             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1152      &         ntermd_1(i,j,k,iblock))
1153 C Martix of D parameters for one dimesional foureir series
1154             do l=1,ntermd_1(i,j,k,iblock)
1155              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1156              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1157              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1158              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1159 c            write(iout,*) "whcodze" ,
1160 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1161             enddo
1162             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1163      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1164      &         v2s(m,l,i,j,k,iblock),
1165      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1166 C Martix of D parameters for two dimesional fourier series
1167             do l=1,ntermd_2(i,j,k,iblock)
1168              do m=1,l-1
1169              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1170              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1171              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1172              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1173              enddo!m
1174             enddo!l
1175           enddo!k
1176         enddo!j
1177       enddo!i
1178       enddo!iblock
1179       if (lprint) then
1180       write (iout,*)
1181       write (iout,*) 'Constants for double torsionals'
1182       do iblock=1,2
1183       do i=0,ntortyp-1
1184         do j=-ntortyp+1,ntortyp-1
1185           do k=-ntortyp+1,ntortyp-1
1186             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1187      &        ' nsingle',ntermd_1(i,j,k,iblock),
1188      &        ' ndouble',ntermd_2(i,j,k,iblock)
1189             write (iout,*)
1190             write (iout,*) 'Single angles:'
1191             do l=1,ntermd_1(i,j,k,iblock)
1192               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1193      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1194      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1195      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1196             enddo
1197             write (iout,*)
1198             write (iout,*) 'Pairs of angles:'
1199             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1200             do l=1,ntermd_2(i,j,k,iblock)
1201               write (iout,'(i5,20f10.5)')
1202      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1203             enddo
1204             write (iout,*)
1205            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1206             do l=1,ntermd_2(i,j,k,iblock)
1207               write (iout,'(i5,20f10.5)')
1208      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1209      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1210             enddo
1211             write (iout,*)
1212           enddo
1213         enddo
1214       enddo
1215       enddo
1216       endif
1217
1218       ELSE IF (TOR_MODE.eq.1) THEN
1219
1220 C read valence-torsional parameters
1221       read (itorp,*,end=113,err=113) ntortyp
1222       nkcctyp=ntortyp
1223       write (iout,*) "Valence-torsional parameters read in ntortyp",
1224      &   ntortyp
1225       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1226       write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1227       do i=-ntyp,-1
1228         itortyp(i)=-itortyp(-i)
1229       enddo
1230       do i=-ntortyp+1,ntortyp-1
1231         do j=-ntortyp+1,ntortyp-1
1232 C first we read the cos and sin gamma parameters
1233           read (itorp,'(13x,a)',end=113,err=113) string
1234           write (iout,*) i,j,string
1235           read (itorp,*,end=113,err=113)
1236      &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1237 C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1238           do k=1,nterm_kcc(j,i)
1239             do l=1,nterm_kcc_Tb(j,i)
1240               do ll=1,nterm_kcc_Tb(j,i)
1241               read (itorp,*,end=113,err=113) ii,jj,kk,
1242      &          v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1243               enddo
1244             enddo
1245           enddo
1246         enddo
1247       enddo
1248
1249       ELSE
1250 c AL 4/8/16: Calculate coefficients from one-body parameters
1251       ntortyp=nloctyp
1252       do i=-ntyp1,ntyp1
1253         itortyp(i)=itype2loc(i)
1254       enddo
1255       write (iout,*)
1256      &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1257      & ,ntortyp
1258       do i=-ntortyp+1,ntortyp-1
1259         do j=-ntortyp+1,ntortyp-1
1260           nterm_kcc(j,i)=2
1261           nterm_kcc_Tb(j,i)=3
1262           do k=1,nterm_kcc_Tb(j,i)
1263             do l=1,nterm_kcc_Tb(j,i)
1264               v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1265      &                         +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1266               v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1267      &                         +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1268             enddo
1269           enddo
1270           do k=1,nterm_kcc_Tb(j,i)
1271             do l=1,nterm_kcc_Tb(j,i)
1272 #ifdef CORRCD
1273               v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1274      &                         -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1275               v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1276      &                         +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1277 #else
1278               v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1279      &                         -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1280               v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1281      &                         +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1282 #endif
1283             enddo
1284           enddo
1285 cf(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1286         enddo
1287       enddo
1288
1289       ENDIF ! TOR_MODE
1290
1291       if (tor_mode.gt.0 .and. lprint) then
1292 c Print valence-torsional parameters
1293         write (iout,'(a)')
1294      &    "Parameters of the valence-torsional potentials"
1295         do i=-ntortyp+1,ntortyp-1
1296         do j=-ntortyp+1,ntortyp-1
1297         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1298         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1299         do k=1,nterm_kcc(j,i)
1300           do l=1,nterm_kcc_Tb(j,i)
1301             do ll=1,nterm_kcc_Tb(j,i)
1302                write (iout,'(3i5,2f15.4)')
1303      &            k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1304             enddo
1305           enddo
1306         enddo
1307         enddo
1308         enddo
1309       endif
1310
1311 #endif
1312 C Read of Side-chain backbone correlation parameters
1313 C Modified 11 May 2012 by Adasko
1314 CCC
1315 C
1316       read (isccor,*,end=119,err=119) nsccortyp
1317       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1318       do i=-ntyp,-1
1319         isccortyp(i)=-isccortyp(-i)
1320       enddo
1321       iscprol=isccortyp(20)
1322 c      write (iout,*) 'ntortyp',ntortyp
1323       maxinter=3
1324 cc maxinter is maximum interaction sites
1325       do l=1,maxinter
1326       do i=1,nsccortyp
1327         do j=1,nsccortyp
1328           read (isccor,*,end=119,err=119)
1329      &nterm_sccor(i,j),nlor_sccor(i,j)
1330 c          write (iout,*) nterm_sccor(i,j)
1331           v0ijsccor=0.0d0
1332           v0ijsccor1=0.0d0
1333           v0ijsccor2=0.0d0
1334           v0ijsccor3=0.0d0
1335           si=-1.0d0
1336           nterm_sccor(-i,j)=nterm_sccor(i,j)
1337           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1338           nterm_sccor(i,-j)=nterm_sccor(i,j)
1339 c          write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1340 c     &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
1341           do k=1,nterm_sccor(i,j)
1342             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1343      &    ,v2sccor(k,l,i,j)
1344             if (j.eq.iscprol) then
1345              if (i.eq.isccortyp(10)) then
1346              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1347              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1348              else
1349              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1350      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1351              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1352      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1353              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1354              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1355              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1356              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1357              endif
1358             else
1359              if (i.eq.isccortyp(10)) then
1360              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1361              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1362              else
1363                if (j.eq.isccortyp(10)) then
1364              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1365              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1366                else
1367              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1368              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1369              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1370              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1371              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1372              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1373                 endif
1374                endif
1375             endif
1376             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1377             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1378             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1379             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1380             si=-si
1381            enddo
1382           do k=1,nlor_sccor(i,j)
1383             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1384      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1385             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1386      &(1+vlor3sccor(k,i,j)**2)
1387           enddo
1388           v0sccor(l,i,j)=v0ijsccor
1389           v0sccor(l,-i,j)=v0ijsccor1
1390           v0sccor(l,i,-j)=v0ijsccor2
1391           v0sccor(l,-i,-j)=v0ijsccor3
1392           enddo
1393         enddo
1394       enddo
1395       close (isccor)
1396       if (lprint) then
1397         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1398         do l=1,maxinter
1399         do i=1,nsccortyp
1400           do j=1,nsccortyp
1401             write (iout,*) 'ityp',i,' jtyp',j
1402             write (iout,*) 'Fourier constants'
1403             do k=1,nterm_sccor(i,j)
1404               write (iout,'(2(1pe15.5))')
1405      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1406             enddo
1407             write (iout,*) 'Lorenz constants'
1408             do k=1,nlor_sccor(i,j)
1409               write (iout,'(3(1pe15.5))')
1410      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1411             enddo
1412           enddo
1413         enddo
1414         enddo
1415       endif
1416 C
1417 C Read electrostatic-interaction parameters
1418 C
1419       if (lprint) then
1420         write (iout,'(/a)') 'Electrostatic interaction constants:'
1421         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1422      &            'IT','JT','APP','BPP','AEL6','AEL3'
1423       endif
1424       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1425       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1426       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1427       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1428       close (ielep)
1429       do i=1,2
1430         do j=1,2
1431         rri=rpp(i,j)**6
1432         app (i,j)=epp(i,j)*rri*rri
1433         bpp (i,j)=-2.0D0*epp(i,j)*rri
1434         ael6(i,j)=elpp6(i,j)*4.2D0**6
1435         ael3(i,j)=elpp3(i,j)*4.2D0**3
1436         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1437      &                    ael6(i,j),ael3(i,j)
1438         enddo
1439       enddo
1440 C
1441 C Read side-chain interaction parameters.
1442 C
1443       read (isidep,*,end=117,err=117) ipot,expon
1444       if (ipot.lt.1 .or. ipot.gt.5) then
1445         write (iout,'(2a)') 'Error while reading SC interaction',
1446      &               'potential file - unknown potential type.'
1447         stop
1448       endif
1449       expon2=expon/2
1450       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1451      & ', exponents are ',expon,2*expon
1452       goto (10,20,30,30,40) ipot
1453 C----------------------- LJ potential ---------------------------------
1454    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1455      &   (sigma0(i),i=1,ntyp)
1456       if (lprint) then
1457         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1458         write (iout,'(a/)') 'The epsilon array:'
1459         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1460         write (iout,'(/a)') 'One-body parameters:'
1461         write (iout,'(a,4x,a)') 'residue','sigma'
1462         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1463       endif
1464       goto 50
1465 C----------------------- LJK potential --------------------------------
1466    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1467      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1468       if (lprint) then
1469         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1470         write (iout,'(a/)') 'The epsilon array:'
1471         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1472         write (iout,'(/a)') 'One-body parameters:'
1473         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1474         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1475      &        i=1,ntyp)
1476       endif
1477       goto 50
1478 C---------------------- GB or BP potential -----------------------------
1479    30 do i=1,ntyp
1480        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1481       enddo
1482       read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1483       read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1484       read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1485       read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1486       do i=1,ntyp
1487        read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1488 C       write(iout,*) "WARNING!!",i,ntyp
1489        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1490 C       do j=1,ntyp
1491 C       epslip(i,j)=epslip(i,j)+0.05d0
1492 C       enddo
1493       enddo
1494 C For the GB potential convert sigma'**2 into chi'
1495       if (ipot.eq.4) then
1496         do i=1,ntyp
1497           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1498         enddo
1499       endif
1500       if (lprint) then
1501         write (iout,'(/a/)') 'Parameters of the BP potential:'
1502         write (iout,'(a/)') 'The epsilon array:'
1503         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1504         write (iout,'(/a)') 'One-body parameters:'
1505         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1506      &       '    chip  ','    alph  '
1507         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1508      &                     chip(i),alp(i),i=1,ntyp)
1509       endif
1510       goto 50
1511 C--------------------- GBV potential -----------------------------------
1512 c   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1513 c     &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1514 c     &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1515    40 do i=1,ntyp
1516        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1517       enddo
1518       read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1519       read (isidep,*,end=117,err=117)(rr0(i),i=1,ntyp)
1520       read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1521       read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1522       read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1523       do i=1,ntyp
1524        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1525        if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1526       enddo
1527       if (lprint) then
1528         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1529         write (iout,'(a/)') 'The epsilon array:'
1530         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1531         write (iout,'(/a)') 'One-body parameters:'
1532         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1533      &      's||/s_|_^2','    chip  ','    alph  '
1534         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1535      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1536       endif
1537    50 continue
1538       close (isidep)
1539 C-----------------------------------------------------------------------
1540 C Calculate the "working" parameters of SC interactions.
1541       do i=2,ntyp
1542         do j=1,i-1
1543           eps(i,j)=eps(j,i)
1544           epslip(i,j)=epslip(j,i)
1545         enddo
1546       enddo
1547       do i=1,ntyp
1548         do j=i,ntyp
1549           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1550           sigma(j,i)=sigma(i,j)
1551           rs0(i,j)=dwa16*sigma(i,j)
1552           rs0(j,i)=rs0(i,j)
1553         enddo
1554       enddo
1555       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1556      & 'Working parameters of the SC interactions:',
1557      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1558      & '  chi1   ','   chi2   '
1559       do i=1,ntyp
1560         do j=i,ntyp
1561           epsij=eps(i,j)
1562           epsijlip=epslip(i,j)
1563           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1564             rrij=sigma(i,j)
1565           else
1566             rrij=rr0(i)+rr0(j)
1567           endif
1568           r0(i,j)=rrij
1569           r0(j,i)=rrij
1570           rrij=rrij**expon
1571           epsij=eps(i,j)
1572           sigeps=dsign(1.0D0,epsij)
1573           epsij=dabs(epsij)
1574           aa_aq(i,j)=epsij*rrij*rrij
1575           bb_aq(i,j)=-sigeps*epsij*rrij
1576           aa_aq(j,i)=aa_aq(i,j)
1577           bb_aq(j,i)=bb_aq(i,j)
1578           sigeps=dsign(1.0D0,epsijlip)
1579           epsijlip=dabs(epsijlip)
1580           aa_lip(i,j)=epsijlip*rrij*rrij
1581           bb_lip(i,j)=-sigeps*epsijlip*rrij
1582           aa_lip(j,i)=aa_lip(i,j)
1583           bb_lip(j,i)=bb_lip(i,j)
1584           if (ipot.gt.2) then
1585             sigt1sq=sigma0(i)**2
1586             sigt2sq=sigma0(j)**2
1587             sigii1=sigii(i)
1588             sigii2=sigii(j)
1589             ratsig1=sigt2sq/sigt1sq
1590             ratsig2=1.0D0/ratsig1
1591             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1592             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1593             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1594           else
1595             rsum_max=sigma(i,j)
1596           endif
1597 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1598             sigmaii(i,j)=rsum_max
1599             sigmaii(j,i)=rsum_max
1600 c         else
1601 c           sigmaii(i,j)=r0(i,j)
1602 c           sigmaii(j,i)=r0(i,j)
1603 c         endif
1604 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1605           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1606             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1607             augm(i,j)=epsij*r_augm**(2*expon)
1608 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1609             augm(j,i)=augm(i,j)
1610           else
1611             augm(i,j)=0.0D0
1612             augm(j,i)=0.0D0
1613           endif
1614           if (lprint) then
1615             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1616      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1617      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1618           endif
1619         enddo
1620       enddo
1621 C
1622 C Define the SC-p interaction constants
1623 C
1624 #ifdef OLDSCP
1625       do i=1,20
1626 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1627 C helix formation)
1628 c       aad(i,1)=0.3D0*4.0D0**12
1629 C Following line for constants currently implemented
1630 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1631         aad(i,1)=1.5D0*4.0D0**12
1632 c       aad(i,1)=0.17D0*5.6D0**12
1633         aad(i,2)=aad(i,1)
1634 C "Soft" SC-p repulsion
1635         bad(i,1)=0.0D0
1636 C Following line for constants currently implemented
1637 c       aad(i,1)=0.3D0*4.0D0**6
1638 C "Hard" SC-p repulsion
1639         bad(i,1)=3.0D0*4.0D0**6
1640 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1641         bad(i,2)=bad(i,1)
1642 c       aad(i,1)=0.0D0
1643 c       aad(i,2)=0.0D0
1644 c       bad(i,1)=1228.8D0
1645 c       bad(i,2)=1228.8D0
1646       enddo
1647 #else
1648 C
1649 C 8/9/01 Read the SC-p interaction constants from file
1650 C
1651       do i=1,ntyp
1652         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1653       enddo
1654       do i=1,ntyp
1655         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1656         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1657         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1658         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1659       enddo
1660
1661       if (lprint) then
1662         write (iout,*) "Parameters of SC-p interactions:"
1663         do i=1,20
1664           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1665      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1666         enddo
1667       endif
1668 #endif
1669 C
1670 C Define the constants of the disulfide bridge
1671 C
1672 C      ebr=-12.0D0
1673 c
1674 c Old arbitrary potential - commented out.
1675 c
1676 c      dbr= 4.20D0
1677 c      fbr= 3.30D0
1678 c
1679 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1680 c energy surface of diethyl disulfide.
1681 c A. Liwo and U. Kozlowska, 11/24/03
1682 c
1683 C      D0CM = 3.78d0
1684 C      AKCM = 15.1d0
1685 C      AKTH = 11.0d0
1686 C      AKCT = 12.0d0
1687 C      V1SS =-1.08d0
1688 C      V2SS = 7.61d0
1689 C      V3SS = 13.7d0
1690       write (iout,*) dyn_ss,'Dynamic disulfide parameters'
1691       if (dyn_ss) then
1692         ss_depth=ebr/wsc-0.25*eps(1,1)
1693         Ht=Ht/wsc-0.25*eps(1,1)
1694         akcm=akcm*wstrain/wsc
1695         akth=akth*wstrain/wsc
1696         akct=akct*wstrain/wsc
1697         v1ss=v1ss*wstrain/wsc
1698         v2ss=v2ss*wstrain/wsc
1699         v3ss=v3ss*wstrain/wsc
1700       else
1701         if (wstrain.ne.0.0) then
1702          ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
1703         else
1704           ss_depth=0.0
1705         endif
1706       endif
1707
1708 C      if (lprint) then
1709       write (iout,'(/a)') "Disulfide bridge parameters:"
1710       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1711       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1712       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1713       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1714      & ' v3ss:',v3ss
1715 C      endif
1716       if (shield_mode.gt.0) then
1717       pi=3.141592d0
1718 C VSolvSphere the volume of solving sphere
1719 C      print *,pi,"pi"
1720 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1721 C there will be no distinction between proline peptide group and normal peptide
1722 C group in case of shielding parameters
1723       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1724       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1725       write (iout,*) "VSolvSphere",VSolvSphere," VSolvSphere_div",
1726      &  VSolvSphere_div
1727 C long axis of side chain
1728       do i=1,ntyp
1729       long_r_sidechain(i)=vbldsc0(1,i)
1730       short_r_sidechain(i)=sigma0(i)
1731       enddo
1732       buff_shield=1.0d0
1733       endif
1734       return
1735   111 write (iout,*) "Error reading bending energy parameters."
1736       goto 999
1737   112 write (iout,*) "Error reading rotamer energy parameters."
1738       goto 999
1739   113 write (iout,*) "Error reading torsional energy parameters."
1740       goto 999
1741   114 write (iout,*) "Error reading double torsional energy parameters."
1742       goto 999
1743   115 write (iout,*)
1744      &  "Error reading cumulant (multibody energy) parameters."
1745       goto 999
1746   116 write (iout,*) "Error reading electrostatic energy parameters."
1747       goto 999
1748  1161 write (iout,*) "Error reading lipid parameters."
1749       goto 999
1750   117 write (iout,*) "Error reading side chain interaction parameters."
1751       goto 999
1752   118 write (iout,*) "Error reading SCp interaction parameters."
1753       goto 999
1754   119 write (iout,*) "Error reading SCCOR parameters"
1755       goto 999
1756   121 write (iout,*) "Error reading bond parameters"
1757   999 continue
1758 #ifdef MPI
1759       call MPI_Finalize(Ierror)
1760 #endif
1761       stop
1762       return
1763       end