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