source/unres/src_MD-M-SAXS/parmread.F

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