source/cluster/wham/src-HCD-5D/parmread.F

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