source/unres/src_MD-M/parmread.F

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