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