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