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