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       do i=-ntyp,-1
 947         isccortyp(i)=-isccortyp(-i)
 948       enddo
 949       maxinter=3
 950 cc maxinter is maximum interaction sites
 951 C NOW READING FOR LL aminoacid chirality
 952       do l=1,maxinter
 953       do i=1,nsccortyp
 954         do j=1,nsccortyp
 955           read (isccor,*,end=119,err=119)
 956      & nterm_sccor(i,j),nlor_sccor(i,j)
 957           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 958           v0ijsccor=0.0d0
 959           si=-1.0d0
 960
 961           do k=1,nterm_sccor(i,j)
 962             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 963      &    ,v2sccor(k,l,i,j)
 964             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 965             si=-si
 966             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 967             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 968
 969           enddo
 970           do k=1,nlor_sccor(i,j)
 971             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 972      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 973             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 974      &(1+vlor3sccor(k,i,j)**2)
 975           enddo
 976           v0sccor(l,i,j)=v0ijsccor
 977           v0sccor(l,-i,-j)=v0ijsccor
 978         enddo
 979       enddo
 980       enddo
 981 C NOW READING FOR LD aminoacid chirality
 982       do l=1,maxinter
 983       do i=1,nsccortyp
 984         do j=-nsccortyp,-1
 985           read (isccor,*,end=119,err=119)
 986      & nterm_sccor(i,j),nlor_sccor(i,j)
 987           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 988           v0ijsccor=0.0d0
 989           si=-1.0d0
 990
 991           do k=1,nterm_sccor(i,j)
 992             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 993      &    ,v2sccor(k,l,i,j)
 994             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 995             si=-si
 996             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 997             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 998
 999           enddo
1000           do k=1,nlor_sccor(i,j)
1001             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1002      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1003             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1004      &(1+vlor3sccor(k,i,j)**2)
1005           enddo
1006           v0sccor(l,i,j)=v0ijsccor
1007           v0sccor(l,-i,-j)=v0ijsccor
1008         enddo
1009       enddo
1010       enddo
1011 C symetry for GLY
1012              i=10
1013              do l=1,maxinter
1014                 do j=-nsccortyp,-1
1015                 do k=1,nterm_sccor(i,j)
1016                 v1sccor(k,l,10,j)=v1sccor(k,l,10,-j)
1017                 v2sccor(k,l,10,j)=-v2sccor(k,l,10,-j)
1018                 v1sccor(k,l,j,10)=v1sccor(k,l,-j,10)
1019                 v2sccor(k,l,j,10)=-v2sccor(k,l,-j,10)
1020                enddo
1021               enddo
1022              enddo
1023
1024       close (isccor)
1025
1026 #endif
1027       if (lprint) then
1028         write (iout,'(/a/)') 'Torsional constants:'
1029         do l=1,maxinter
1030         do i=1,nsccortyp
1031           do j=1,nsccortyp
1032             write (iout,*) 'ityp',i,' jtyp',j
1033             write (iout,*) 'Fourier constants'
1034             do k=1,nterm_sccor(i,j)
1035       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1036             enddo
1037             write (iout,*) 'Lorenz constants'
1038             do k=1,nlor_sccor(i,j)
1039               write (iout,'(3(1pe15.5))')
1040      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1041             enddo
1042           enddo
1043         enddo
1044         enddo
1045       endif
1046
1047 C
1048 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1049 C         interaction energy of the Gly, Ala, and Pro prototypes.
1050 C
1051       if (lprint) then
1052         write (iout,*)
1053         write (iout,*) "Coefficients of the cumulants"
1054       endif
1055       read (ifourier,*) nloctyp
1056 #ifdef NEWCORR
1057       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1058       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1059       itype2loc(ntyp1)=nloctyp
1060       iloctyp(nloctyp)=ntyp1
1061 #else
1062       do i=1,ntyp1
1063         itype2loc(i)=itortyp(i)
1064       enddo
1065       iloctyp(0)=10
1066       iloctyp(1)=9
1067       iloctyp(2)=20
1068       iloctyp(3)=ntyp1
1069 #endif
1070       do i=1,ntyp1
1071         itype2loc(-i)=-itype2loc(i)
1072       enddo
1073       do i=1,nloctyp
1074         iloctyp(-i)=-iloctyp(i)
1075       enddo
1076       write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1077       write (iout,*) "nloctyp",nloctyp,
1078      &  " iloctyp",(iloctyp(i),i=0,nloctyp)
1079       do i=0,nloctyp-1
1080         read (ifourier,*,end=115,err=115)
1081         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1082 #ifdef NEWCORR
1083         read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
1084         read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
1085         read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
1086         read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
1087         read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
1088 #endif
1089         if (lprint) then
1090         write (iout,*) 'Type',i
1091         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1092         endif
1093         b(3,-i)= b(3,i)
1094         b(5,-i)=-b(5,i)
1095         b(4,-i)=-b(4,i)
1096         b(2,-i)= b(2,i)
1097         B1(1,i)  = b(3,i)
1098         B1(2,i)  = b(5,i)
1099         B1(1,-i) = b(3,i)
1100         B1(2,-i) = -b(5,i)
1101 c        b1(1,i)=0.0d0
1102 c        b1(2,i)=0.0d0
1103         B1tilde(1,i) = b(3,i)
1104         B1tilde(2,i) =-b(5,i)
1105         B1tilde(1,-i) =-b(3,i)
1106         B1tilde(2,-i) =b(5,i)
1107 c        b1tilde(1,i)=0.0d0
1108 c        b1tilde(2,i)=0.0d0
1109         B2(1,i)  = b(2,i)
1110         B2(2,i)  = b(4,i)
1111         B2(1,-i)  =b(2,i)
1112         B2(2,-i)  =-b(4,i)
1113 cc        B1tilde(1,i) = b(3,i)
1114 cc        B1tilde(2,i) =-b(5,i)
1115 C        B1tilde(1,-i) =-b(3,i)
1116 C        B1tilde(2,-i) =b(5,i)
1117 cc        b1tilde(1,i)=0.0d0
1118 cc        b1tilde(2,i)=0.0d0
1119 cc        B2(1,i)  = b(2,i)
1120 cc        B2(2,i)  = b(4,i)
1121 C        B2(1,-i)  =b(2,i)
1122 C        B2(2,-i)  =-b(4,i)
1123
1124 c        b2(1,i)=0.0d0
1125 c        b2(2,i)=0.0d0
1126         CC(1,1,i)= b(7,i)
1127         CC(2,2,i)=-b(7,i)
1128         CC(2,1,i)= b(9,i)
1129         CC(1,2,i)= b(9,i)
1130         CC(1,1,-i)= b(7,i)
1131         CC(2,2,-i)=-b(7,i)
1132         CC(2,1,-i)=-b(9,i)
1133         CC(1,2,-i)=-b(9,i)
1134 c        CC(1,1,i)=0.0d0
1135 c        CC(2,2,i)=0.0d0
1136 c        CC(2,1,i)=0.0d0
1137 c        CC(1,2,i)=0.0d0
1138         Ctilde(1,1,i)=b(7,i)
1139         Ctilde(1,2,i)=b(9,i)
1140         Ctilde(2,1,i)=-b(9,i)
1141         Ctilde(2,2,i)=b(7,i)
1142         Ctilde(1,1,-i)=b(7,i)
1143         Ctilde(1,2,-i)=-b(9,i)
1144         Ctilde(2,1,-i)=b(9,i)
1145         Ctilde(2,2,-i)=b(7,i)
1146
1147 c        Ctilde(1,1,i)=0.0d0
1148 c        Ctilde(1,2,i)=0.0d0
1149 c        Ctilde(2,1,i)=0.0d0
1150 c        Ctilde(2,2,i)=0.0d0
1151         DD(1,1,i)= b(6,i)
1152         DD(2,2,i)=-b(6,i)
1153         DD(2,1,i)= b(8,i)
1154         DD(1,2,i)= b(8,i)
1155         DD(1,1,-i)= b(6,i)
1156         DD(2,2,-i)=-b(6,i)
1157         DD(2,1,-i)=-b(8,i)
1158         DD(1,2,-i)=-b(8,i)
1159 c        DD(1,1,i)=0.0d0
1160 c        DD(2,2,i)=0.0d0
1161 c        DD(2,1,i)=0.0d0
1162 c        DD(1,2,i)=0.0d0
1163         Dtilde(1,1,i)=b(6,i)
1164         Dtilde(1,2,i)=b(8,i)
1165         Dtilde(2,1,i)=-b(8,i)
1166         Dtilde(2,2,i)=b(6,i)
1167         Dtilde(1,1,-i)=b(6,i)
1168         Dtilde(1,2,-i)=-b(8,i)
1169         Dtilde(2,1,-i)=b(8,i)
1170         Dtilde(2,2,-i)=b(6,i)
1171
1172 c        Dtilde(1,1,i)=0.0d0
1173 c        Dtilde(1,2,i)=0.0d0
1174 c        Dtilde(2,1,i)=0.0d0
1175 c        Dtilde(2,2,i)=0.0d0
1176         EEold(1,1,i)= b(10,i)+b(11,i)
1177         EEold(2,2,i)=-b(10,i)+b(11,i)
1178         EEold(2,1,i)= b(12,i)-b(13,i)
1179         EEold(1,2,i)= b(12,i)+b(13,i)
1180         EEold(1,1,-i)= b(10,i)+b(11,i)
1181         EEold(2,2,-i)=-b(10,i)+b(11,i)
1182         EEold(2,1,-i)=-b(12,i)+b(13,i)
1183         EEold(1,2,-i)=-b(12,i)-b(13,i)
1184         write(iout,*) "TU DOCHODZE"
1185         print *,"JESTEM"
1186 c        ee(1,1,i)=1.0d0
1187 c        ee(2,2,i)=1.0d0
1188 c        ee(2,1,i)=0.0d0
1189 c        ee(1,2,i)=0.0d0
1190 c        ee(2,1,i)=ee(1,2,i)
1191       enddo
1192       lprint=.true.
1193       if (lprint) then
1194       do i=1,nloctyp
1195         write (iout,*) 'Type',i
1196         write (iout,*) 'B1'
1197         write(iout,*) B1(1,i),B1(2,i)
1198         write (iout,*) 'B2'
1199         write(iout,*) B2(1,i),B2(2,i)
1200         write (iout,*) 'CC'
1201         do j=1,2
1202           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1203         enddo
1204         write(iout,*) 'DD'
1205         do j=1,2
1206           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1207         enddo
1208         write(iout,*) 'EE'
1209         do j=1,2
1210           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1211         enddo
1212       enddo
1213       endif
1214       lprint=.false.
1215
1216 C
1217 C Read electrostatic-interaction parameters
1218 C
1219       if (lprint) then
1220         write (iout,*)
1221         write (iout,'(/a)') 'Electrostatic interaction constants:'
1222         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1223      &            'IT','JT','APP','BPP','AEL6','AEL3'
1224       endif
1225       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1226       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1227       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1228       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1229       close (ielep)
1230       do i=1,2
1231         do j=1,2
1232         rri=rpp(i,j)**6
1233         app (i,j)=epp(i,j)*rri*rri
1234         bpp (i,j)=-2.0D0*epp(i,j)*rri
1235         ael6(i,j)=elpp6(i,j)*4.2D0**6
1236         ael3(i,j)=elpp3(i,j)*4.2D0**3
1237 c        lprint=.true.
1238         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1239      &                    ael6(i,j),ael3(i,j)
1240 c        lprint=.false.
1241         enddo
1242       enddo
1243 C
1244 C Read side-chain interaction parameters.
1245 C
1246       read (isidep,*,end=117,err=117) ipot,expon
1247       if (ipot.lt.1 .or. ipot.gt.5) then
1248         write (iout,'(2a)') 'Error while reading SC interaction',
1249      &               'potential file - unknown potential type.'
1250 #ifdef MPI
1251         call MPI_Finalize(Ierror)
1252 #endif
1253         stop
1254       endif
1255       expon2=expon/2
1256       if(me.eq.king)
1257      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1258      & ', exponents are ',expon,2*expon
1259       goto (10,20,30,30,40) ipot
1260 C----------------------- LJ potential ---------------------------------
1261    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1262      &   (sigma0(i),i=1,ntyp)
1263       if (lprint) then
1264         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1265         write (iout,'(a/)') 'The epsilon array:'
1266         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1267         write (iout,'(/a)') 'One-body parameters:'
1268         write (iout,'(a,4x,a)') 'residue','sigma'
1269         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1270       endif
1271       goto 50
1272 C----------------------- LJK potential --------------------------------
1273    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1274      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1275       if (lprint) then
1276         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1277         write (iout,'(a/)') 'The epsilon array:'
1278         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1279         write (iout,'(/a)') 'One-body parameters:'
1280         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1281         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1282      &        i=1,ntyp)
1283       endif
1284       goto 50
1285 C---------------------- GB or BP potential -----------------------------
1286    30 do i=1,ntyp
1287        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1288       enddo
1289       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1290       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1291       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1292       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1293 C now we start reading lipid
1294       do i=1,ntyp
1295        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1296        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1297
1298 C       print *,"WARNING!!"
1299 C       do j=1,ntyp
1300 C       epslip(i,j)=epslip(i,j)+0.05d0
1301 C       enddo
1302       enddo
1303       write(iout,*) epslip(1,1),"OK?"
1304 C For the GB potential convert sigma'**2 into chi'
1305       if (ipot.eq.4) then
1306         do i=1,ntyp
1307           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1308         enddo
1309       endif
1310       if (lprint) then
1311         write (iout,'(/a/)') 'Parameters of the BP potential:'
1312         write (iout,'(a/)') 'The epsilon array:'
1313         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1314         write (iout,'(/a)') 'One-body parameters:'
1315         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1316      &       '    chip  ','    alph  '
1317         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1318      &                     chip(i),alp(i),i=1,ntyp)
1319       endif
1320       goto 50
1321 C--------------------- GBV potential -----------------------------------
1322    40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1323      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1324      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1325       if (lprint) then
1326         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1327         write (iout,'(a/)') 'The epsilon array:'
1328         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1329         write (iout,'(/a)') 'One-body parameters:'
1330         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1331      &      's||/s_|_^2','    chip  ','    alph  '
1332         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1333      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1334       endif
1335    50 continue
1336       close (isidep)
1337 C-----------------------------------------------------------------------
1338 C Calculate the "working" parameters of SC interactions.
1339       do i=2,ntyp
1340         do j=1,i-1
1341           eps(i,j)=eps(j,i)
1342           epslip(i,j)=epslip(j,i)
1343         enddo
1344       enddo
1345       do i=1,ntyp
1346         do j=i,ntyp
1347           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1348           sigma(j,i)=sigma(i,j)
1349           rs0(i,j)=dwa16*sigma(i,j)
1350           rs0(j,i)=rs0(i,j)
1351         enddo
1352       enddo
1353       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1354      & 'Working parameters of the SC interactions:',
1355      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1356      & '  chi1   ','   chi2   '
1357       do i=1,ntyp
1358         do j=i,ntyp
1359           epsij=eps(i,j)
1360           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1361             rrij=sigma(i,j)
1362           else
1363             rrij=rr0(i)+rr0(j)
1364           endif
1365           r0(i,j)=rrij
1366           r0(j,i)=rrij
1367           rrij=rrij**expon
1368           epsij=eps(i,j)
1369           sigeps=dsign(1.0D0,epsij)
1370           epsij=dabs(epsij)
1371           aa_aq(i,j)=epsij*rrij*rrij
1372           bb_aq(i,j)=-sigeps*epsij*rrij
1373           aa_aq(j,i)=aa_aq(i,j)
1374           bb_aq(j,i)=bb_aq(i,j)
1375           epsijlip=epslip(i,j)
1376           sigeps=dsign(1.0D0,epsijlip)
1377           epsijlip=dabs(epsijlip)
1378           aa_lip(i,j)=epsijlip*rrij*rrij
1379           bb_lip(i,j)=-sigeps*epsijlip*rrij
1380           aa_lip(j,i)=aa_lip(i,j)
1381           bb_lip(j,i)=bb_lip(i,j)
1382           if (ipot.gt.2) then
1383             sigt1sq=sigma0(i)**2
1384             sigt2sq=sigma0(j)**2
1385             sigii1=sigii(i)
1386             sigii2=sigii(j)
1387             ratsig1=sigt2sq/sigt1sq
1388             ratsig2=1.0D0/ratsig1
1389             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1390             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1391             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1392           else
1393             rsum_max=sigma(i,j)
1394           endif
1395 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1396             sigmaii(i,j)=rsum_max
1397             sigmaii(j,i)=rsum_max
1398 c         else
1399 c           sigmaii(i,j)=r0(i,j)
1400 c           sigmaii(j,i)=r0(i,j)
1401 c         endif
1402 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1403           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1404             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1405             augm(i,j)=epsij*r_augm**(2*expon)
1406 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1407             augm(j,i)=augm(i,j)
1408           else
1409             augm(i,j)=0.0D0
1410             augm(j,i)=0.0D0
1411           endif
1412           if (lprint) then
1413             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1414      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1415      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1416           endif
1417         enddo
1418       enddo
1419       write(iout,*) "tube param"
1420       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep,
1421      & ccavtubpep,dcavtubpep,tubetranenepep
1422       sigmapeptube=sigmapeptube**6
1423       sigeps=dsign(1.0D0,epspeptube)
1424       epspeptube=dabs(epspeptube)
1425       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1426       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1427       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1428       do i=1,ntyp
1429        read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i),
1430      & ccavtub(i),dcavtub(i),tubetranene(i)
1431        sigmasctube=sigmasctube**6
1432        sigeps=dsign(1.0D0,epssctube)
1433        epssctube=dabs(epssctube)
1434        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1435        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1436       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1437       enddo
1438
1439 #ifdef OLDSCP
1440 C
1441 C Define the SC-p interaction constants (hard-coded; old style)
1442 C
1443       do i=1,ntyp
1444 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1445 C helix formation)
1446 c       aad(i,1)=0.3D0*4.0D0**12
1447 C Following line for constants currently implemented
1448 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1449         aad(i,1)=1.5D0*4.0D0**12
1450 c       aad(i,1)=0.17D0*5.6D0**12
1451         aad(i,2)=aad(i,1)
1452 C "Soft" SC-p repulsion
1453         bad(i,1)=0.0D0
1454 C Following line for constants currently implemented
1455 c       aad(i,1)=0.3D0*4.0D0**6
1456 C "Hard" SC-p repulsion
1457         bad(i,1)=3.0D0*4.0D0**6
1458 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1459         bad(i,2)=bad(i,1)
1460 c       aad(i,1)=0.0D0
1461 c       aad(i,2)=0.0D0
1462 c       bad(i,1)=1228.8D0
1463 c       bad(i,2)=1228.8D0
1464       enddo
1465 #else
1466 C
1467 C 8/9/01 Read the SC-p interaction constants from file
1468 C
1469       do i=1,ntyp
1470         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1471       enddo
1472       do i=1,ntyp
1473         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1474         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1475         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1476         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1477       enddo
1478 c      lprint=.true.
1479       if (lprint) then
1480         write (iout,*) "Parameters of SC-p interactions:"
1481         do i=1,ntyp
1482           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1483      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1484         enddo
1485       endif
1486 c      lprint=.false.
1487 #endif
1488 C
1489 C Define the constants of the disulfide bridge
1490 C
1491 C      ebr=-12.00D0
1492 c
1493 c Old arbitrary potential - commented out.
1494 c
1495 c      dbr= 4.20D0
1496 c      fbr= 3.30D0
1497 c
1498 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1499 c energy surface of diethyl disulfide.
1500 c A. Liwo and U. Kozlowska, 11/24/03
1501 c
1502 C      D0CM = 3.78d0
1503 C      AKCM = 15.1d0
1504 C      AKTH = 11.0d0
1505 C      AKCT = 12.0d0
1506 C      V1SS =-1.08d0
1507 C      V2SS = 7.61d0
1508 C      V3SS = 13.7d0
1509 c      akcm=0.0d0
1510 c      akth=0.0d0
1511 c      akct=0.0d0
1512 c      v1ss=0.0d0
1513 c      v2ss=0.0d0
1514 c      v3ss=0.0d0
1515
1516 C      if(me.eq.king) then
1517 C      write (iout,'(/a)') "Disulfide bridge parameters:"
1518 C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1519 C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1520 C      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1521 C      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1522 C     &  ' v3ss:',v3ss
1523 C      endif
1524 C set the variables used for shielding effect
1525 C      write (iout,*) "SHIELD MODE",shield_mode
1526 C      if (shield_mode.gt.0) then
1527 C VSolvSphere the volume of solving sphere
1528 C      print *,pi,"pi"
1529 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1530 C there will be no distinction between proline peptide group and normal peptide
1531 C group in case of shielding parameters
1532 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1533 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1534 C      write (iout,*) VSolvSphere,VSolvSphere_div
1535 C long axis of side chain
1536 C      do i=1,ntyp
1537 C      long_r_sidechain(i)=vbldsc0(1,i)
1538 C      short_r_sidechain(i)=sigma0(i)
1539 C      enddo
1540 C lets set the buffor value
1541 C      buff_shield=1.0d0
1542 C      endif
1543       return
1544   111 write (iout,*) "Error reading bending energy parameters."
1545       goto 999
1546   112 write (iout,*) "Error reading rotamer energy parameters."
1547       goto 999
1548   113 write (iout,*) "Error reading torsional energy parameters."
1549       goto 999
1550   114 write (iout,*) "Error reading double torsional energy parameters."
1551       goto 999
1552   115 write (iout,*)
1553      &  "Error reading cumulant (multibody energy) parameters."
1554       goto 999
1555   116 write (iout,*) "Error reading electrostatic energy parameters."
1556       goto 999
1557  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1558       goto 999
1559   117 write (iout,*) "Error reading side chain interaction parameters."
1560       goto 999
1561   118 write (iout,*) "Error reading SCp interaction parameters."
1562       goto 999
1563   119 write (iout,*) "Error reading SCCOR parameters"
1564       go to 999
1565   121 write (iout,*) "Error in Czybyshev parameters"
1566   999 continue
1567 #ifdef MPI
1568       call MPI_Finalize(Ierror)
1569 #endif
1570       stop
1571       return
1572       end
1573
1574
1575       subroutine getenv_loc(var, val)
1576       character(*) var, val
1577
1578 #ifdef WINIFL
1579       character(2000) line
1580       external ilen
1581
1582       open (196,file='env',status='old',readonly,shared)
1583       iread=0
1584 c      write(*,*)'looking for ',var
1585 10    read(196,*,err=11,end=11)line
1586       iread=index(line,var)
1587 c      write(*,*)iread,' ',var,' ',line
1588       if (iread.eq.0) go to 10
1589 c      write(*,*)'---> ',line
1590 11    continue
1591       if(iread.eq.0) then
1592 c       write(*,*)'CHUJ'
1593        val=''
1594       else
1595        iread=iread+ilen(var)+1
1596        read (line(iread:),*,err=12,end=12) val
1597 c       write(*,*)'OK: ',var,' = ',val
1598       endif
1599       close(196)
1600       return
1601 12    val=''
1602       close(196)
1603 #elif (defined CRAY)
1604       integer lennam,lenval,ierror
1605 c
1606 c        getenv using a POSIX call, useful on the T3D
1607 c        Sept 1996, comment out error check on advice of H. Pritchard
1608 c
1609       lennam = len(var)
1610       if(lennam.le.0) stop '--error calling getenv--'
1611       call pxfgetenv(var,lennam,val,lenval,ierror)
1612 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1613 #else
1614       call getenv(var,val)
1615 #endif
1616 C set the variables used for shielding effect
1617 C      if (shield_mode.gt.0) then
1618 C VSolvSphere the volume of solving sphere
1619 C      print *,pi,"pi"
1620 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1621 C there will be no distinction between proline peptide group and normal peptide
1622 C group in case of shielding parameters
1623 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1624 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1625 C long axis of side chain
1626 C      do i=1,ntyp
1627 C      long_r_sidechain(i)=vbldsc0(1,i)
1628 C      short_r_sidechain(i)=sigma0(i)
1629 C      enddo
1630 C lets set the buffor value
1631 C      buff_shield=1.0d0
1632 C      endif
1633       return
1634       end