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       itortyp(ntyp1)=ntortyp
 634       do iblock=1,2
 635       do i=-ntyp,-1
 636        itortyp(i)=-itortyp(-i)
 637       enddo
 638       write (iout,*) 'ntortyp',ntortyp
 639       do i=0,ntortyp-1
 640         do j=-ntortyp+1,ntortyp-1
 641           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 642      &          nlor(i,j,iblock)
 643           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 644           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 645           v0ij=0.0d0
 646           si=-1.0d0
 647           do k=1,nterm(i,j,iblock)
 648             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 649      &      v2(k,i,j,iblock)
 650             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 651             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 652             v0ij=v0ij+si*v1(k,i,j,iblock)
 653             si=-si
 654 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
 655 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
 656 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
 657           enddo
 658           do k=1,nlor(i,j,iblock)
 659             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 660      &        vlor2(k,i,j),vlor3(k,i,j)
 661             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 662           enddo
 663           v0(i,j,iblock)=v0ij
 664           v0(-i,-j,iblock)=v0ij
 665         enddo
 666       enddo
 667       enddo
 668       close (itorp)
 669       if (lprint) then
 670         write (iout,'(/a/)') 'Torsional constants:'
 671         do iblock=1,2
 672         do i=0,ntortyp-1
 673           do j=-ntortyp+1,ntortyp-1
 674             write (iout,*) 'ityp',i,' jtyp',j
 675             write (iout,*) 'Fourier constants'
 676             do k=1,nterm(i,j,iblock)
 677               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 678      &        v2(k,i,j,iblock)
 679             enddo
 680             write (iout,*) 'Lorenz constants'
 681             do k=1,nlor(i,j,iblock)
 682               write (iout,'(3(1pe15.5))')
 683      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 684             enddo
 685           enddo
 686         enddo
 687         enddo
 688       endif
 689
 690 C
 691 C 6/23/01 Read parameters for double torsionals
 692 C
 693       do iblock=1,2
 694       do i=0,ntortyp-1
 695         do j=-ntortyp+1,ntortyp-1
 696           do k=-ntortyp+1,ntortyp-1
 697             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 698 c              write (iout,*) "OK onelett",
 699 c     &         i,j,k,t1,t2,t3
 700
 701             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 702      &        .or. t3.ne.toronelet(k)) then
 703               write (iout,*) "Error in double torsional parameter file",
 704      &         i,j,k,t1,t2,t3
 705 #ifdef MPI
 706               call MPI_Finalize(Ierror)
 707 #endif
 708                stop "Error in double torsional parameter file"
 709             endif
 710            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
 711      &         ntermd_2(i,j,k,iblock)
 712             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 713             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 714             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
 715      &         ntermd_1(i,j,k,iblock))
 716             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
 717      &         ntermd_1(i,j,k,iblock))
 718             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
 719      &         ntermd_1(i,j,k,iblock))
 720             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
 721      &         ntermd_1(i,j,k,iblock))
 722 C Martix of D parameters for one dimesional foureir series
 723             do l=1,ntermd_1(i,j,k,iblock)
 724              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 725              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 726              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 727              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 728 c            write(iout,*) "whcodze" ,
 729 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 730             enddo
 731             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
 732      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 733      &         v2s(m,l,i,j,k,iblock),
 734      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 735 C Martix of D parameters for two dimesional fourier series
 736             do l=1,ntermd_2(i,j,k,iblock)
 737              do m=1,l-1
 738              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 739              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 740              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 741              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 742              enddo!m
 743             enddo!l
 744           enddo!k
 745         enddo!j
 746       enddo!i
 747       enddo!iblock
 748       if (lprint) then
 749       write (iout,*)
 750       write (iout,*) 'Constants for double torsionals'
 751       do iblock=1,2
 752       do i=0,ntortyp-1
 753         do j=-ntortyp+1,ntortyp-1
 754           do k=-ntortyp+1,ntortyp-1
 755             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 756      &        ' nsingle',ntermd_1(i,j,k,iblock),
 757      &        ' ndouble',ntermd_2(i,j,k,iblock)
 758             write (iout,*)
 759             write (iout,*) 'Single angles:'
 760             do l=1,ntermd_1(i,j,k,iblock)
 761               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 762      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 763      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 764      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 765             enddo
 766             write (iout,*)
 767             write (iout,*) 'Pairs of angles:'
 768             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 769             do l=1,ntermd_2(i,j,k,iblock)
 770               write (iout,'(i5,20f10.5)')
 771      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 772             enddo
 773             write (iout,*)
 774             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 775             do l=1,ntermd_2(i,j,k,iblock)
 776               write (iout,'(i5,20f10.5)')
 777      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 778      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 779             enddo
 780             write (iout,*)
 781           enddo
 782         enddo
 783       enddo
 784       enddo
 785       endif
 786 #endif
 787 C read Czybyshev torsional parameters
 788       read (itorkcc,*,end=121,err=121) nkcctyp
 789       read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
 790       do i=-ntyp,-1
 791         itortyp_kcc(i)=-itortyp_kcc(-i)
 792       enddo
 793       do i=0,nkcctyp
 794         do j=0,nkcctyp
 795 C first we read the cos and sin gamma parameters
 796           read (itorkcc,*,end=121,err=121)
 797      &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
 798 C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
 799           do k=1,nterm_kcc(j,i)
 800             do l=1,nterm_kcc_Tb(j,i)
 801               read (itorkcc,*,end=121,err=121) v11_chyb(l,k,j,i)
 802             enddo
 803             do l=1,nterm_kcc_Tb(j,i)
 804               read (itorkcc,*,end=121,err=121) v21_chyb(l,k,j,i)
 805             enddo
 806             do l=1,nterm_kcc_Tb(j,i)
 807               read (itorkcc,*,end=121,err=121) v12_chyb(l,k,j,i)
 808             enddo
 809             do l=1,nterm_kcc_Tb(j,i)
 810               read (itorkcc,*,end=121,err=121) v22_chyb(l,k,j,i)
 811             enddo
 812             read (itorkcc,*,end=121,err=121) v1_kcc(k,j,i)
 813             read (itorkcc,*,end=121,err=121) v2_kcc(k,j,i)
 814           enddo
 815         enddo
 816       enddo
 817       if (lprint) then
 818 c Print valence-torsional parameters
 819         write (iout,'(a)')
 820      &    "Parameters of the valence-torsional potentials"
 821         do i=0,nkcctyp
 822         do j=0,nkcctyp
 823         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
 824         write (iout,'(2a20,a15)') "v_kcc","v1_chyb","v2_chyb"
 825         do k=1,nterm_kcc(j,i)
 826           write (iout,'(i5,f15.10,i5,2f15.10)')
 827      &      k,v1_kcc(k,j,i),1,v11_chyb(1,k,j,i),v21_chyb(1,k,j,i)
 828           do l=2,nterm_kcc_Tb(j,i)
 829             write (iout,'(20x,i5,2f15.10)')
 830      &        l,v11_chyb(l,k,j,i),v21_chyb(l,k,j,i)
 831           enddo
 832           write (iout,'(i5,f15.10,i5,2f15.10)')
 833      &      k,v2_kcc(k,j,i),1,v12_chyb(1,k,j,i),v22_chyb(1,k,j,i)
 834           do l=2,nterm_kcc_Tb(j,i)
 835             write (iout,'(20x,i5,2f15.10)')
 836      &        l,v12_chyb(l,k,j,i),v22_chyb(l,k,j,i)
 837           enddo
 838           write (iout,'(a)')
 839         enddo
 840         enddo
 841         enddo
 842       endif
 843 C here will be the apropriate recalibrating for D-aminoacid
 844 C        read (ithetkcc,*,end=121,err=121) nkcctyp
 845       do i=0,nkcctyp
 846         read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i)
 847         do j=1,nbend_kcc_Tb(i)
 848           read (ithetkcc,*,end=121,err=121) v1bend_chyb(j,i)
 849         enddo
 850       enddo
 851       if (lprint) then
 852         write (iout,'(a)')
 853      &    "Parameters of the valence-only potentials"
 854         do i=0,nkcctyp
 855         write (iout,'(2a)') "Type ",toronelet(i)
 856         do k=1,nbend_kcc_Tb(i)
 857           write(iout,'(i5,f15.10)') k,v1bend_chyb(k,i)
 858         enddo
 859         enddo
 860       endif
 861 C Read of Side-chain backbone correlation parameters
 862 C Modified 11 May 2012 by Adasko
 863 CCC
 864 C
 865       read (isccor,*,end=119,err=119) nsccortyp
 866 #ifdef SCCORPDB
 867       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 868       do i=-ntyp,-1
 869         isccortyp(i)=-isccortyp(-i)
 870       enddo
 871       iscprol=isccortyp(20)
 872 c      write (iout,*) 'ntortyp',ntortyp
 873       maxinter=3
 874 cc maxinter is maximum interaction sites
 875       do l=1,maxinter
 876       do i=1,nsccortyp
 877         do j=1,nsccortyp
 878           read (isccor,*,end=119,err=119)
 879      &nterm_sccor(i,j),nlor_sccor(i,j)
 880           v0ijsccor=0.0d0
 881           v0ijsccor1=0.0d0
 882           v0ijsccor2=0.0d0
 883           v0ijsccor3=0.0d0
 884           si=-1.0d0
 885           nterm_sccor(-i,j)=nterm_sccor(i,j)
 886           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 887           nterm_sccor(i,-j)=nterm_sccor(i,j)
 888           do k=1,nterm_sccor(i,j)
 889             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 890      &    ,v2sccor(k,l,i,j)
 891             if (j.eq.iscprol) then
 892              if (i.eq.isccortyp(10)) then
 893              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 894              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 895              else
 896              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 897      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 898              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 899      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 900              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 901              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 902              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 903              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 904              endif
 905             else
 906              if (i.eq.isccortyp(10)) then
 907              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 908              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 909              else
 910                if (j.eq.isccortyp(10)) then
 911              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 912              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 913                else
 914              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 915              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 916              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 917              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 918              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 919              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 920                 endif
 921                endif
 922             endif
 923             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 924             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 925             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 926             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 927             si=-si
 928           enddo
 929           do k=1,nlor_sccor(i,j)
 930             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 931      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 932             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 933      &(1+vlor3sccor(k,i,j)**2)
 934           enddo
 935           v0sccor(l,i,j)=v0ijsccor
 936           v0sccor(l,-i,j)=v0ijsccor1
 937           v0sccor(l,i,-j)=v0ijsccor2
 938           v0sccor(l,-i,-j)=v0ijsccor3
 939         enddo
 940       enddo
 941       enddo
 942       close (isccor)
 943 #else
 944       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 945 c      write (iout,*) 'ntortyp',ntortyp
 946       maxinter=3
 947 cc maxinter is maximum interaction sites
 948       do l=1,maxinter
 949       do i=1,nsccortyp
 950         do j=1,nsccortyp
 951           read (isccor,*,end=119,err=119)
 952      & nterm_sccor(i,j),nlor_sccor(i,j)
 953           v0ijsccor=0.0d0
 954           si=-1.0d0
 955
 956           do k=1,nterm_sccor(i,j)
 957             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 958      &    ,v2sccor(k,l,i,j)
 959             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 960             si=-si
 961           enddo
 962           do k=1,nlor_sccor(i,j)
 963             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 964      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 965             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 966      &(1+vlor3sccor(k,i,j)**2)
 967           enddo
 968           v0sccor(l,i,j)=v0ijsccor
 969         enddo
 970       enddo
 971       enddo
 972       close (isccor)
 973
 974 #endif
 975       if (lprint) then
 976         write (iout,'(/a/)') 'Torsional constants:'
 977         do l=1,maxinter
 978         do i=1,nsccortyp
 979           do j=1,nsccortyp
 980             write (iout,*) 'ityp',i,' jtyp',j
 981             write (iout,*) 'Fourier constants'
 982             do k=1,nterm_sccor(i,j)
 983       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 984             enddo
 985             write (iout,*) 'Lorenz constants'
 986             do k=1,nlor_sccor(i,j)
 987               write (iout,'(3(1pe15.5))')
 988      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 989             enddo
 990           enddo
 991         enddo
 992         enddo
 993       endif
 994
 995 C
 996 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 997 C         interaction energy of the Gly, Ala, and Pro prototypes.
 998 C
 999       if (lprint) then
1000         write (iout,*)
1001         write (iout,*) "Coefficients of the cumulants"
1002       endif
1003       read (ifourier,*) nloctyp
1004 #ifdef NEWCORR
1005       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1006       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1007       itype2loc(ntyp1)=nloctyp
1008       iloctyp(nloctyp)=ntyp1
1009 #else
1010       do i=1,ntyp1
1011         itype2loc(i)=itortyp(i)
1012       enddo
1013       iloctyp(0)=10
1014       iloctyp(1)=9
1015       iloctyp(2)=20
1016       iloctyp(3)=ntyp1
1017 #endif
1018       do i=1,ntyp1
1019         itype2loc(-i)=-itype2loc(i)
1020       enddo
1021       do i=1,nloctyp
1022         iloctyp(-i)=-iloctyp(i)
1023       enddo
1024       write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1025       write (iout,*) "nloctyp",nloctyp,
1026      &  " iloctyp",(iloctyp(i),i=0,nloctyp)
1027       do i=0,nloctyp-1
1028         read (ifourier,*,end=115,err=115)
1029         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1030 #ifdef NEWCORR
1031         read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
1032         read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
1033         read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
1034         read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
1035         read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
1036 #endif
1037         if (lprint) then
1038         write (iout,*) 'Type',i
1039         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1040         endif
1041         b(3,-i)= b(3,i)
1042         b(5,-i)=-b(5,i)
1043         b(4,-i)=-b(4,i)
1044         b(2,-i)= b(2,i)
1045         B1(1,i)  = b(3,i)
1046         B1(2,i)  = b(5,i)
1047         B1(1,-i) = b(3,i)
1048         B1(2,-i) = -b(5,i)
1049 c        b1(1,i)=0.0d0
1050 c        b1(2,i)=0.0d0
1051         B1tilde(1,i) = b(3,i)
1052         B1tilde(2,i) =-b(5,i)
1053         B1tilde(1,-i) =-b(3,i)
1054         B1tilde(2,-i) =b(5,i)
1055 c        b1tilde(1,i)=0.0d0
1056 c        b1tilde(2,i)=0.0d0
1057         B2(1,i)  = b(2,i)
1058         B2(2,i)  = b(4,i)
1059         B2(1,-i)  =b(2,i)
1060         B2(2,-i)  =-b(4,i)
1061 cc        B1tilde(1,i) = b(3,i)
1062 cc        B1tilde(2,i) =-b(5,i)
1063 C        B1tilde(1,-i) =-b(3,i)
1064 C        B1tilde(2,-i) =b(5,i)
1065 cc        b1tilde(1,i)=0.0d0
1066 cc        b1tilde(2,i)=0.0d0
1067 cc        B2(1,i)  = b(2,i)
1068 cc        B2(2,i)  = b(4,i)
1069 C        B2(1,-i)  =b(2,i)
1070 C        B2(2,-i)  =-b(4,i)
1071
1072 c        b2(1,i)=0.0d0
1073 c        b2(2,i)=0.0d0
1074         CC(1,1,i)= b(7,i)
1075         CC(2,2,i)=-b(7,i)
1076         CC(2,1,i)= b(9,i)
1077         CC(1,2,i)= b(9,i)
1078         CC(1,1,-i)= b(7,i)
1079         CC(2,2,-i)=-b(7,i)
1080         CC(2,1,-i)=-b(9,i)
1081         CC(1,2,-i)=-b(9,i)
1082 c        CC(1,1,i)=0.0d0
1083 c        CC(2,2,i)=0.0d0
1084 c        CC(2,1,i)=0.0d0
1085 c        CC(1,2,i)=0.0d0
1086         Ctilde(1,1,i)=b(7,i)
1087         Ctilde(1,2,i)=b(9,i)
1088         Ctilde(2,1,i)=-b(9,i)
1089         Ctilde(2,2,i)=b(7,i)
1090         Ctilde(1,1,-i)=b(7,i)
1091         Ctilde(1,2,-i)=-b(9,i)
1092         Ctilde(2,1,-i)=b(9,i)
1093         Ctilde(2,2,-i)=b(7,i)
1094
1095 c        Ctilde(1,1,i)=0.0d0
1096 c        Ctilde(1,2,i)=0.0d0
1097 c        Ctilde(2,1,i)=0.0d0
1098 c        Ctilde(2,2,i)=0.0d0
1099         DD(1,1,i)= b(6,i)
1100         DD(2,2,i)=-b(6,i)
1101         DD(2,1,i)= b(8,i)
1102         DD(1,2,i)= b(8,i)
1103         DD(1,1,-i)= b(6,i)
1104         DD(2,2,-i)=-b(6,i)
1105         DD(2,1,-i)=-b(8,i)
1106         DD(1,2,-i)=-b(8,i)
1107 c        DD(1,1,i)=0.0d0
1108 c        DD(2,2,i)=0.0d0
1109 c        DD(2,1,i)=0.0d0
1110 c        DD(1,2,i)=0.0d0
1111         Dtilde(1,1,i)=b(6,i)
1112         Dtilde(1,2,i)=b(8,i)
1113         Dtilde(2,1,i)=-b(8,i)
1114         Dtilde(2,2,i)=b(6,i)
1115         Dtilde(1,1,-i)=b(6,i)
1116         Dtilde(1,2,-i)=-b(8,i)
1117         Dtilde(2,1,-i)=b(8,i)
1118         Dtilde(2,2,-i)=b(6,i)
1119
1120 c        Dtilde(1,1,i)=0.0d0
1121 c        Dtilde(1,2,i)=0.0d0
1122 c        Dtilde(2,1,i)=0.0d0
1123 c        Dtilde(2,2,i)=0.0d0
1124         EEold(1,1,i)= b(10,i)+b(11,i)
1125         EEold(2,2,i)=-b(10,i)+b(11,i)
1126         EEold(2,1,i)= b(12,i)-b(13,i)
1127         EEold(1,2,i)= b(12,i)+b(13,i)
1128         EEold(1,1,-i)= b(10,i)+b(11,i)
1129         EEold(2,2,-i)=-b(10,i)+b(11,i)
1130         EEold(2,1,-i)=-b(12,i)+b(13,i)
1131         EEold(1,2,-i)=-b(12,i)-b(13,i)
1132         write(iout,*) "TU DOCHODZE"
1133         print *,"JESTEM"
1134 c        ee(1,1,i)=1.0d0
1135 c        ee(2,2,i)=1.0d0
1136 c        ee(2,1,i)=0.0d0
1137 c        ee(1,2,i)=0.0d0
1138 c        ee(2,1,i)=ee(1,2,i)
1139       enddo
1140       lprint=.true.
1141       if (lprint) then
1142       do i=1,nloctyp
1143         write (iout,*) 'Type',i
1144         write (iout,*) 'B1'
1145         write(iout,*) B1(1,i),B1(2,i)
1146         write (iout,*) 'B2'
1147         write(iout,*) B2(1,i),B2(2,i)
1148         write (iout,*) 'CC'
1149         do j=1,2
1150           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1151         enddo
1152         write(iout,*) 'DD'
1153         do j=1,2
1154           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1155         enddo
1156         write(iout,*) 'EE'
1157         do j=1,2
1158           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1159         enddo
1160       enddo
1161       endif
1162       lprint=.false.
1163
1164 C
1165 C Read electrostatic-interaction parameters
1166 C
1167       if (lprint) then
1168         write (iout,*)
1169         write (iout,'(/a)') 'Electrostatic interaction constants:'
1170         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1171      &            'IT','JT','APP','BPP','AEL6','AEL3'
1172       endif
1173       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1174       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1175       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1176       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1177       close (ielep)
1178       do i=1,2
1179         do j=1,2
1180         rri=rpp(i,j)**6
1181         app (i,j)=epp(i,j)*rri*rri
1182         bpp (i,j)=-2.0D0*epp(i,j)*rri
1183         ael6(i,j)=elpp6(i,j)*4.2D0**6
1184         ael3(i,j)=elpp3(i,j)*4.2D0**3
1185 c        lprint=.true.
1186         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1187      &                    ael6(i,j),ael3(i,j)
1188 c        lprint=.false.
1189         enddo
1190       enddo
1191 C
1192 C Read side-chain interaction parameters.
1193 C
1194       read (isidep,*,end=117,err=117) ipot,expon
1195       if (ipot.lt.1 .or. ipot.gt.5) then
1196         write (iout,'(2a)') 'Error while reading SC interaction',
1197      &               'potential file - unknown potential type.'
1198 #ifdef MPI
1199         call MPI_Finalize(Ierror)
1200 #endif
1201         stop
1202       endif
1203       expon2=expon/2
1204       if(me.eq.king)
1205      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1206      & ', exponents are ',expon,2*expon
1207       goto (10,20,30,30,40) ipot
1208 C----------------------- LJ potential ---------------------------------
1209    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1210      &   (sigma0(i),i=1,ntyp)
1211       if (lprint) then
1212         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1213         write (iout,'(a/)') 'The epsilon array:'
1214         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1215         write (iout,'(/a)') 'One-body parameters:'
1216         write (iout,'(a,4x,a)') 'residue','sigma'
1217         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1218       endif
1219       goto 50
1220 C----------------------- LJK potential --------------------------------
1221    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1222      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1223       if (lprint) then
1224         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1225         write (iout,'(a/)') 'The epsilon array:'
1226         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1227         write (iout,'(/a)') 'One-body parameters:'
1228         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1229         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1230      &        i=1,ntyp)
1231       endif
1232       goto 50
1233 C---------------------- GB or BP potential -----------------------------
1234    30 do i=1,ntyp
1235        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1236       enddo
1237       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1238       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1239       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1240       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1241 C now we start reading lipid
1242       do i=1,ntyp
1243        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1244        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1245
1246 C       print *,"WARNING!!"
1247 C       do j=1,ntyp
1248 C       epslip(i,j)=epslip(i,j)+0.05d0
1249 C       enddo
1250       enddo
1251       write(iout,*) epslip(1,1),"OK?"
1252 C For the GB potential convert sigma'**2 into chi'
1253       if (ipot.eq.4) then
1254         do i=1,ntyp
1255           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1256         enddo
1257       endif
1258       if (lprint) then
1259         write (iout,'(/a/)') 'Parameters of the BP potential:'
1260         write (iout,'(a/)') 'The epsilon array:'
1261         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1262         write (iout,'(/a)') 'One-body parameters:'
1263         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1264      &       '    chip  ','    alph  '
1265         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1266      &                     chip(i),alp(i),i=1,ntyp)
1267       endif
1268       goto 50
1269 C--------------------- GBV potential -----------------------------------
1270    40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1271      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1272      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1273       if (lprint) then
1274         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1275         write (iout,'(a/)') 'The epsilon array:'
1276         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1277         write (iout,'(/a)') 'One-body parameters:'
1278         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1279      &      's||/s_|_^2','    chip  ','    alph  '
1280         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1281      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1282       endif
1283    50 continue
1284       close (isidep)
1285 C-----------------------------------------------------------------------
1286 C Calculate the "working" parameters of SC interactions.
1287       do i=2,ntyp
1288         do j=1,i-1
1289           eps(i,j)=eps(j,i)
1290           epslip(i,j)=epslip(j,i)
1291         enddo
1292       enddo
1293       do i=1,ntyp
1294         do j=i,ntyp
1295           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1296           sigma(j,i)=sigma(i,j)
1297           rs0(i,j)=dwa16*sigma(i,j)
1298           rs0(j,i)=rs0(i,j)
1299         enddo
1300       enddo
1301       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1302      & 'Working parameters of the SC interactions:',
1303      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1304      & '  chi1   ','   chi2   '
1305       do i=1,ntyp
1306         do j=i,ntyp
1307           epsij=eps(i,j)
1308           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1309             rrij=sigma(i,j)
1310           else
1311             rrij=rr0(i)+rr0(j)
1312           endif
1313           r0(i,j)=rrij
1314           r0(j,i)=rrij
1315           rrij=rrij**expon
1316           epsij=eps(i,j)
1317           sigeps=dsign(1.0D0,epsij)
1318           epsij=dabs(epsij)
1319           aa_aq(i,j)=epsij*rrij*rrij
1320           bb_aq(i,j)=-sigeps*epsij*rrij
1321           aa_aq(j,i)=aa_aq(i,j)
1322           bb_aq(j,i)=bb_aq(i,j)
1323           epsijlip=epslip(i,j)
1324           sigeps=dsign(1.0D0,epsijlip)
1325           epsijlip=dabs(epsijlip)
1326           aa_lip(i,j)=epsijlip*rrij*rrij
1327           bb_lip(i,j)=-sigeps*epsijlip*rrij
1328           aa_lip(j,i)=aa_lip(i,j)
1329           bb_lip(j,i)=bb_lip(i,j)
1330           if (ipot.gt.2) then
1331             sigt1sq=sigma0(i)**2
1332             sigt2sq=sigma0(j)**2
1333             sigii1=sigii(i)
1334             sigii2=sigii(j)
1335             ratsig1=sigt2sq/sigt1sq
1336             ratsig2=1.0D0/ratsig1
1337             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1338             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1339             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1340           else
1341             rsum_max=sigma(i,j)
1342           endif
1343 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1344             sigmaii(i,j)=rsum_max
1345             sigmaii(j,i)=rsum_max
1346 c         else
1347 c           sigmaii(i,j)=r0(i,j)
1348 c           sigmaii(j,i)=r0(i,j)
1349 c         endif
1350 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1351           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1352             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1353             augm(i,j)=epsij*r_augm**(2*expon)
1354 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1355             augm(j,i)=augm(i,j)
1356           else
1357             augm(i,j)=0.0D0
1358             augm(j,i)=0.0D0
1359           endif
1360           if (lprint) then
1361             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1362      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1363      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1364           endif
1365         enddo
1366       enddo
1367       write(iout,*) "tube param"
1368       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep,
1369      & ccavtubpep,dcavtubpep,tubetranenepep
1370       sigmapeptube=sigmapeptube**6
1371       sigeps=dsign(1.0D0,epspeptube)
1372       epspeptube=dabs(epspeptube)
1373       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1374       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1375       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1376       do i=1,ntyp
1377        read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i),
1378      & ccavtub(i),dcavtub(i),tubetranene(i)
1379        sigmasctube=sigmasctube**6
1380        sigeps=dsign(1.0D0,epssctube)
1381        epssctube=dabs(epssctube)
1382        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1383        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1384       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1385       enddo
1386
1387 #ifdef OLDSCP
1388 C
1389 C Define the SC-p interaction constants (hard-coded; old style)
1390 C
1391       do i=1,ntyp
1392 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1393 C helix formation)
1394 c       aad(i,1)=0.3D0*4.0D0**12
1395 C Following line for constants currently implemented
1396 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1397         aad(i,1)=1.5D0*4.0D0**12
1398 c       aad(i,1)=0.17D0*5.6D0**12
1399         aad(i,2)=aad(i,1)
1400 C "Soft" SC-p repulsion
1401         bad(i,1)=0.0D0
1402 C Following line for constants currently implemented
1403 c       aad(i,1)=0.3D0*4.0D0**6
1404 C "Hard" SC-p repulsion
1405         bad(i,1)=3.0D0*4.0D0**6
1406 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1407         bad(i,2)=bad(i,1)
1408 c       aad(i,1)=0.0D0
1409 c       aad(i,2)=0.0D0
1410 c       bad(i,1)=1228.8D0
1411 c       bad(i,2)=1228.8D0
1412       enddo
1413 #else
1414 C
1415 C 8/9/01 Read the SC-p interaction constants from file
1416 C
1417       do i=1,ntyp
1418         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1419       enddo
1420       do i=1,ntyp
1421         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1422         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1423         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1424         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1425       enddo
1426 c      lprint=.true.
1427       if (lprint) then
1428         write (iout,*) "Parameters of SC-p interactions:"
1429         do i=1,ntyp
1430           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1431      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1432         enddo
1433       endif
1434 c      lprint=.false.
1435 #endif
1436 C
1437 C Define the constants of the disulfide bridge
1438 C
1439 C      ebr=-12.00D0
1440 c
1441 c Old arbitrary potential - commented out.
1442 c
1443 c      dbr= 4.20D0
1444 c      fbr= 3.30D0
1445 c
1446 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1447 c energy surface of diethyl disulfide.
1448 c A. Liwo and U. Kozlowska, 11/24/03
1449 c
1450 C      D0CM = 3.78d0
1451 C      AKCM = 15.1d0
1452 C      AKTH = 11.0d0
1453 C      AKCT = 12.0d0
1454 C      V1SS =-1.08d0
1455 C      V2SS = 7.61d0
1456 C      V3SS = 13.7d0
1457 c      akcm=0.0d0
1458 c      akth=0.0d0
1459 c      akct=0.0d0
1460 c      v1ss=0.0d0
1461 c      v2ss=0.0d0
1462 c      v3ss=0.0d0
1463
1464 C      if(me.eq.king) then
1465 C      write (iout,'(/a)') "Disulfide bridge parameters:"
1466 C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1467 C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1468 C      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1469 C      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1470 C     &  ' v3ss:',v3ss
1471 C      endif
1472 C set the variables used for shielding effect
1473 C      write (iout,*) "SHIELD MODE",shield_mode
1474 C      if (shield_mode.gt.0) then
1475 C VSolvSphere the volume of solving sphere
1476 C      print *,pi,"pi"
1477 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1478 C there will be no distinction between proline peptide group and normal peptide
1479 C group in case of shielding parameters
1480 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1481 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1482 C      write (iout,*) VSolvSphere,VSolvSphere_div
1483 C long axis of side chain
1484 C      do i=1,ntyp
1485 C      long_r_sidechain(i)=vbldsc0(1,i)
1486 C      short_r_sidechain(i)=sigma0(i)
1487 C      enddo
1488 C lets set the buffor value
1489 C      buff_shield=1.0d0
1490 C      endif
1491       return
1492   111 write (iout,*) "Error reading bending energy parameters."
1493       goto 999
1494   112 write (iout,*) "Error reading rotamer energy parameters."
1495       goto 999
1496   113 write (iout,*) "Error reading torsional energy parameters."
1497       goto 999
1498   114 write (iout,*) "Error reading double torsional energy parameters."
1499       goto 999
1500   115 write (iout,*)
1501      &  "Error reading cumulant (multibody energy) parameters."
1502       goto 999
1503   116 write (iout,*) "Error reading electrostatic energy parameters."
1504       goto 999
1505  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1506       goto 999
1507   117 write (iout,*) "Error reading side chain interaction parameters."
1508       goto 999
1509   118 write (iout,*) "Error reading SCp interaction parameters."
1510       goto 999
1511   119 write (iout,*) "Error reading SCCOR parameters"
1512       go to 999
1513   121 write (iout,*) "Error in Czybyshev parameters"
1514   999 continue
1515 #ifdef MPI
1516       call MPI_Finalize(Ierror)
1517 #endif
1518       stop
1519       return
1520       end
1521
1522
1523       subroutine getenv_loc(var, val)
1524       character(*) var, val
1525
1526 #ifdef WINIFL
1527       character(2000) line
1528       external ilen
1529
1530       open (196,file='env',status='old',readonly,shared)
1531       iread=0
1532 c      write(*,*)'looking for ',var
1533 10    read(196,*,err=11,end=11)line
1534       iread=index(line,var)
1535 c      write(*,*)iread,' ',var,' ',line
1536       if (iread.eq.0) go to 10
1537 c      write(*,*)'---> ',line
1538 11    continue
1539       if(iread.eq.0) then
1540 c       write(*,*)'CHUJ'
1541        val=''
1542       else
1543        iread=iread+ilen(var)+1
1544        read (line(iread:),*,err=12,end=12) val
1545 c       write(*,*)'OK: ',var,' = ',val
1546       endif
1547       close(196)
1548       return
1549 12    val=''
1550       close(196)
1551 #elif (defined CRAY)
1552       integer lennam,lenval,ierror
1553 c
1554 c        getenv using a POSIX call, useful on the T3D
1555 c        Sept 1996, comment out error check on advice of H. Pritchard
1556 c
1557       lennam = len(var)
1558       if(lennam.le.0) stop '--error calling getenv--'
1559       call pxfgetenv(var,lennam,val,lenval,ierror)
1560 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1561 #else
1562       call getenv(var,val)
1563 #endif
1564 C set the variables used for shielding effect
1565 C      if (shield_mode.gt.0) then
1566 C VSolvSphere the volume of solving sphere
1567 C      print *,pi,"pi"
1568 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1569 C there will be no distinction between proline peptide group and normal peptide
1570 C group in case of shielding parameters
1571 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1572 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1573 C long axis of side chain
1574 C      do i=1,ntyp
1575 C      long_r_sidechain(i)=vbldsc0(1,i)
1576 C      short_r_sidechain(i)=sigma0(i)
1577 C      enddo
1578 C lets set the buffor value
1579 C      buff_shield=1.0d0
1580 C      endif
1581       return
1582       end