source/unres/src_MD-M-SAXS.safe/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    40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1554      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1555      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1556       if (lprint) then
1557         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1558         write (iout,'(a/)') 'The epsilon array:'
1559         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1560         write (iout,'(/a)') 'One-body parameters:'
1561         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1562      &      's||/s_|_^2','    chip  ','    alph  '
1563         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1564      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1565       endif
1566    50 continue
1567       close (isidep)
1568 C-----------------------------------------------------------------------
1569 C Calculate the "working" parameters of SC interactions.
1570       do i=2,ntyp
1571         do j=1,i-1
1572           eps(i,j)=eps(j,i)
1573           epslip(i,j)=epslip(j,i)
1574         enddo
1575       enddo
1576       do i=1,ntyp
1577         do j=i,ntyp
1578           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1579           sigma(j,i)=sigma(i,j)
1580           rs0(i,j)=dwa16*sigma(i,j)
1581           rs0(j,i)=rs0(i,j)
1582         enddo
1583       enddo
1584       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1585      & 'Working parameters of the SC interactions:',
1586      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1587      & '  chi1   ','   chi2   '
1588       do i=1,ntyp
1589         do j=i,ntyp
1590           epsij=eps(i,j)
1591           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1592             rrij=sigma(i,j)
1593           else
1594             rrij=rr0(i)+rr0(j)
1595           endif
1596           r0(i,j)=rrij
1597           r0(j,i)=rrij
1598           rrij=rrij**expon
1599           epsij=eps(i,j)
1600           sigeps=dsign(1.0D0,epsij)
1601           epsij=dabs(epsij)
1602           aa_aq(i,j)=epsij*rrij*rrij
1603           bb_aq(i,j)=-sigeps*epsij*rrij
1604           aa_aq(j,i)=aa_aq(i,j)
1605           bb_aq(j,i)=bb_aq(i,j)
1606           epsijlip=epslip(i,j)
1607           sigeps=dsign(1.0D0,epsijlip)
1608           epsijlip=dabs(epsijlip)
1609           aa_lip(i,j)=epsijlip*rrij*rrij
1610           bb_lip(i,j)=-sigeps*epsijlip*rrij
1611           aa_lip(j,i)=aa_lip(i,j)
1612           bb_lip(j,i)=bb_lip(i,j)
1613           if (ipot.gt.2) then
1614             sigt1sq=sigma0(i)**2
1615             sigt2sq=sigma0(j)**2
1616             sigii1=sigii(i)
1617             sigii2=sigii(j)
1618             ratsig1=sigt2sq/sigt1sq
1619             ratsig2=1.0D0/ratsig1
1620             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1621             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1622             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1623           else
1624             rsum_max=sigma(i,j)
1625           endif
1626 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1627             sigmaii(i,j)=rsum_max
1628             sigmaii(j,i)=rsum_max
1629 c         else
1630 c           sigmaii(i,j)=r0(i,j)
1631 c           sigmaii(j,i)=r0(i,j)
1632 c         endif
1633 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1634           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1635             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1636             augm(i,j)=epsij*r_augm**(2*expon)
1637 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1638             augm(j,i)=augm(i,j)
1639           else
1640             augm(i,j)=0.0D0
1641             augm(j,i)=0.0D0
1642           endif
1643           if (lprint) then
1644             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1645      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1646      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1647           endif
1648         enddo
1649       enddo
1650 #ifdef TUBE
1651       write(iout,*) "tube param"
1652       read(itube,*) epspeptube,sigmapeptube
1653       sigmapeptube=sigmapeptube**6
1654       sigeps=dsign(1.0D0,epspeptube)
1655       epspeptube=dabs(epspeptube)
1656       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1657       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1658       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1659       do i=1,ntyp
1660        read(itube,*) epssctube,sigmasctube
1661        sigmasctube=sigmasctube**6
1662        sigeps=dsign(1.0D0,epssctube)
1663        epssctube=dabs(epssctube)
1664        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1665        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1666       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1667       enddo
1668 #endif
1669
1670 #ifdef OLDSCP
1671 C
1672 C Define the SC-p interaction constants (hard-coded; old style)
1673 C
1674       do i=1,ntyp
1675 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1676 C helix formation)
1677 c       aad(i,1)=0.3D0*4.0D0**12
1678 C Following line for constants currently implemented
1679 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1680         aad(i,1)=1.5D0*4.0D0**12
1681 c       aad(i,1)=0.17D0*5.6D0**12
1682         aad(i,2)=aad(i,1)
1683 C "Soft" SC-p repulsion
1684         bad(i,1)=0.0D0
1685 C Following line for constants currently implemented
1686 c       aad(i,1)=0.3D0*4.0D0**6
1687 C "Hard" SC-p repulsion
1688         bad(i,1)=3.0D0*4.0D0**6
1689 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1690         bad(i,2)=bad(i,1)
1691 c       aad(i,1)=0.0D0
1692 c       aad(i,2)=0.0D0
1693 c       bad(i,1)=1228.8D0
1694 c       bad(i,2)=1228.8D0
1695       enddo
1696 #else
1697 C
1698 C 8/9/01 Read the SC-p interaction constants from file
1699 C
1700       do i=1,ntyp
1701         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1702       enddo
1703       do i=1,ntyp
1704         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1705         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1706         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1707         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1708       enddo
1709 c      lprint=.true.
1710       if (lprint) then
1711         write (iout,*) "Parameters of SC-p interactions:"
1712         do i=1,ntyp
1713           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1714      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1715         enddo
1716       endif
1717 c      lprint=.false.
1718 #endif
1719 C
1720 C Define the constants of the disulfide bridge
1721 C
1722 C      ebr=-12.00D0
1723 c
1724 c Old arbitrary potential - commented out.
1725 c
1726 c      dbr= 4.20D0
1727 c      fbr= 3.30D0
1728 c
1729 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1730 c energy surface of diethyl disulfide.
1731 c A. Liwo and U. Kozlowska, 11/24/03
1732 c
1733 C      D0CM = 3.78d0
1734 C      AKCM = 15.1d0
1735 C      AKTH = 11.0d0
1736 C      AKCT = 12.0d0
1737 C      V1SS =-1.08d0
1738 C      V2SS = 7.61d0
1739 C      V3SS = 13.7d0
1740 c      akcm=0.0d0
1741 c      akth=0.0d0
1742 c      akct=0.0d0
1743 c      v1ss=0.0d0
1744 c      v2ss=0.0d0
1745 c      v3ss=0.0d0
1746
1747 C      if(me.eq.king) then
1748 C      write (iout,'(/a)') "Disulfide bridge parameters:"
1749 C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1750 C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1751 C      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1752 C      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1753 C     &  ' v3ss:',v3ss
1754 C      endif
1755       if (shield_mode.gt.0) then
1756 C VSolvSphere the volume of solving sphere
1757 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1758 C there will be no distinction between proline peptide group and normal peptide
1759 C group in case of shielding parameters
1760 c      write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
1761       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1762       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1763       write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
1764      &  VSolvSphere_div
1765 C long axis of side chain
1766       do i=1,ntyp
1767       long_r_sidechain(i)=vbldsc0(1,i)
1768       short_r_sidechain(i)=sigma0(i)
1769       enddo
1770       buff_shield=1.0d0
1771       endif
1772       return
1773   111 write (iout,*) "Error reading bending energy parameters."
1774       goto 999
1775   112 write (iout,*) "Error reading rotamer energy parameters."
1776       goto 999
1777   113 write (iout,*) "Error reading torsional energy parameters."
1778       goto 999
1779   114 write (iout,*) "Error reading double torsional energy parameters."
1780       goto 999
1781   115 write (iout,*)
1782      &  "Error reading cumulant (multibody energy) parameters."
1783       goto 999
1784   116 write (iout,*) "Error reading electrostatic energy parameters."
1785       goto 999
1786  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1787       goto 999
1788   117 write (iout,*) "Error reading side chain interaction parameters."
1789       goto 999
1790   118 write (iout,*) "Error reading SCp interaction parameters."
1791       goto 999
1792   119 write (iout,*) "Error reading SCCOR parameters"
1793       go to 999
1794   121 write (iout,*) "Error in Czybyshev parameters"
1795   999 continue
1796 #ifdef MPI
1797       call MPI_Finalize(Ierror)
1798 #endif
1799       stop
1800       return
1801       end
1802
1803
1804       subroutine getenv_loc(var, val)
1805       character(*) var, val
1806
1807 #ifdef WINIFL
1808       character(2000) line
1809       external ilen
1810
1811       open (196,file='env',status='old',readonly,shared)
1812       iread=0
1813 c      write(*,*)'looking for ',var
1814 10    read(196,*,err=11,end=11)line
1815       iread=index(line,var)
1816 c      write(*,*)iread,' ',var,' ',line
1817       if (iread.eq.0) go to 10
1818 c      write(*,*)'---> ',line
1819 11    continue
1820       if(iread.eq.0) then
1821 c       write(*,*)'CHUJ'
1822        val=''
1823       else
1824        iread=iread+ilen(var)+1
1825        read (line(iread:),*,err=12,end=12) val
1826 c       write(*,*)'OK: ',var,' = ',val
1827       endif
1828       close(196)
1829       return
1830 12    val=''
1831       close(196)
1832 #elif (defined CRAY)
1833       integer lennam,lenval,ierror
1834 c
1835 c        getenv using a POSIX call, useful on the T3D
1836 c        Sept 1996, comment out error check on advice of H. Pritchard
1837 c
1838       lennam = len(var)
1839       if(lennam.le.0) stop '--error calling getenv--'
1840       call pxfgetenv(var,lennam,val,lenval,ierror)
1841 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1842 #else
1843       call getenv(var,val)
1844 #endif
1845 C set the variables used for shielding effect
1846 C      if (shield_mode.gt.0) then
1847 C VSolvSphere the volume of solving sphere
1848 C      print *,pi,"pi"
1849 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1850 C there will be no distinction between proline peptide group and normal peptide
1851 C group in case of shielding parameters
1852 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1853 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1854 C long axis of side chain
1855 C      do i=1,ntyp
1856 C      long_r_sidechain(i)=vbldsc0(1,i)
1857 C      short_r_sidechain(i)=sigma0(i)
1858 C      enddo
1859 C lets set the buffor value
1860 C      buff_shield=1.0d0
1861 C      endif
1862       return
1863       end