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