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

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