source/unres/src-HCD-5D/parmread.F

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