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