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