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       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),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),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)
 951         CC(2,2,i)=-b(7)
 952         CC(2,1,i)= b(9)
 953         CC(1,2,i)= b(9)
 954         CC(1,1,-i)= b(7)
 955         CC(2,2,-i)=-b(7)
 956         CC(2,1,-i)=-b(9)
 957         CC(1,2,-i)=-b(9)
 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)
 963         Ctilde(1,2,i)=b(9)
 964         Ctilde(2,1,i)=-b(9)
 965         Ctilde(2,2,i)=b(7)
 966         Ctilde(1,1,-i)=b(7)
 967         Ctilde(1,2,-i)=-b(9)
 968         Ctilde(2,1,-i)=b(9)
 969         Ctilde(2,2,-i)=b(7)
 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)
 976         DD(2,2,i)=-b(6)
 977         DD(2,1,i)= b(8)
 978         DD(1,2,i)= b(8)
 979         DD(1,1,-i)= b(6)
 980         DD(2,2,-i)=-b(6)
 981         DD(2,1,-i)=-b(8)
 982         DD(1,2,-i)=-b(8)
 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)
 988         Dtilde(1,2,i)=b(8)
 989         Dtilde(2,1,i)=-b(8)
 990         Dtilde(2,2,i)=b(6)
 991         Dtilde(1,1,-i)=b(6)
 992         Dtilde(1,2,-i)=-b(8)
 993         Dtilde(2,1,-i)=b(8)
 994         Dtilde(2,2,-i)=b(6)
 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)+b(11)
1001         EEold(2,2,i)=-b(10)+b(11)
1002         EEold(2,1,i)= b(12)-b(13)
1003         EEold(1,2,i)= b(12)+b(13)
1004         EEold(1,1,-i)= b(10)+b(11)
1005         EEold(2,2,-i)=-b(10)+b(11)
1006         EEold(2,1,-i)=-b(12)+b(13)
1007         EEold(1,2,-i)=-b(12)-b(13)
1008         write(iout,*) "TU DOCHODZE"
1009 c        ee(1,1,i)=1.0d0
1010 c        ee(2,2,i)=1.0d0
1011 c        ee(2,1,i)=0.0d0
1012 c        ee(1,2,i)=0.0d0
1013 c        ee(2,1,i)=ee(1,2,i)
1014       enddo
1015 c      lprint=.true.
1016       if (lprint) then
1017       do i=1,nloctyp
1018         write (iout,*) 'Type',i
1019         write (iout,*) 'B1'
1020         write(iout,*) B1(1,i),B1(2,i)
1021         write (iout,*) 'B2'
1022         write(iout,*) B2(1,i),B2(2,i)
1023         write (iout,*) 'CC'
1024         do j=1,2
1025           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1026         enddo
1027         write(iout,*) 'DD'
1028         do j=1,2
1029           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1030         enddo
1031         write(iout,*) 'EE'
1032         do j=1,2
1033           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1034         enddo
1035       enddo
1036       endif
1037 c      lprint=.false.
1038
1039 C
1040 C Read electrostatic-interaction parameters
1041 C
1042       if (lprint) then
1043         write (iout,*)
1044         write (iout,'(/a)') 'Electrostatic interaction constants:'
1045         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1046      &            'IT','JT','APP','BPP','AEL6','AEL3'
1047       endif
1048       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1049       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1050       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1051       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1052       close (ielep)
1053       do i=1,2
1054         do j=1,2
1055         rri=rpp(i,j)**6
1056         app (i,j)=epp(i,j)*rri*rri
1057         bpp (i,j)=-2.0D0*epp(i,j)*rri
1058         ael6(i,j)=elpp6(i,j)*4.2D0**6
1059         ael3(i,j)=elpp3(i,j)*4.2D0**3
1060 c        lprint=.true.
1061         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1062      &                    ael6(i,j),ael3(i,j)
1063 c        lprint=.false.
1064         enddo
1065       enddo
1066 C
1067 C Read side-chain interaction parameters.
1068 C
1069       read (isidep,*,end=117,err=117) ipot,expon
1070       if (ipot.lt.1 .or. ipot.gt.5) then
1071         write (iout,'(2a)') 'Error while reading SC interaction',
1072      &               'potential file - unknown potential type.'
1073 #ifdef MPI
1074         call MPI_Finalize(Ierror)
1075 #endif
1076         stop
1077       endif
1078       expon2=expon/2
1079       if(me.eq.king)
1080      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1081      & ', exponents are ',expon,2*expon
1082       goto (10,20,30,30,40) ipot
1083 C----------------------- LJ potential ---------------------------------
1084    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1085      &   (sigma0(i),i=1,ntyp)
1086       if (lprint) then
1087         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1088         write (iout,'(a/)') 'The epsilon array:'
1089         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1090         write (iout,'(/a)') 'One-body parameters:'
1091         write (iout,'(a,4x,a)') 'residue','sigma'
1092         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1093       endif
1094       goto 50
1095 C----------------------- LJK potential --------------------------------
1096    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1097      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1098       if (lprint) then
1099         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1100         write (iout,'(a/)') 'The epsilon array:'
1101         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1102         write (iout,'(/a)') 'One-body parameters:'
1103         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1104         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1105      &        i=1,ntyp)
1106       endif
1107       goto 50
1108 C---------------------- GB or BP potential -----------------------------
1109    30 do i=1,ntyp
1110        read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1111       enddo
1112       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1113       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1114       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1115       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1116 C now we start reading lipid
1117       do i=1,ntyp
1118        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1119        print *,"WARNING!!"
1120        do j=1,ntyp
1121        epslip(i,j)=epslip(i,j)+0.5d0
1122        enddo
1123       enddo
1124 C For the GB potential convert sigma'**2 into chi'
1125       if (ipot.eq.4) then
1126         do i=1,ntyp
1127           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1128         enddo
1129       endif
1130       if (lprint) then
1131         write (iout,'(/a/)') 'Parameters of the BP potential:'
1132         write (iout,'(a/)') 'The epsilon array:'
1133         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1134         write (iout,'(/a)') 'One-body parameters:'
1135         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1136      &       '    chip  ','    alph  '
1137         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1138      &                     chip(i),alp(i),i=1,ntyp)
1139       endif
1140       goto 50
1141 C--------------------- GBV potential -----------------------------------
1142    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1143      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1144      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1145       if (lprint) then
1146         write (iout,'(/a/)') 'Parameters of the GBV 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,5a)') 'residue','   sigma  ','    r0    ',
1151      &      's||/s_|_^2','    chip  ','    alph  '
1152         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1153      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1154       endif
1155    50 continue
1156       close (isidep)
1157 C-----------------------------------------------------------------------
1158 C Calculate the "working" parameters of SC interactions.
1159       do i=2,ntyp
1160         do j=1,i-1
1161           eps(i,j)=eps(j,i)
1162         enddo
1163       enddo
1164       do i=1,ntyp
1165         do j=i,ntyp
1166           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1167           sigma(j,i)=sigma(i,j)
1168           rs0(i,j)=dwa16*sigma(i,j)
1169           rs0(j,i)=rs0(i,j)
1170         enddo
1171       enddo
1172       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1173      & 'Working parameters of the SC interactions:',
1174      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1175      & '  chi1   ','   chi2   '
1176       do i=1,ntyp
1177         do j=i,ntyp
1178           epsij=eps(i,j)
1179           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1180             rrij=sigma(i,j)
1181           else
1182             rrij=rr0(i)+rr0(j)
1183           endif
1184           r0(i,j)=rrij
1185           r0(j,i)=rrij
1186           rrij=rrij**expon
1187           epsij=eps(i,j)
1188           sigeps=dsign(1.0D0,epsij)
1189           epsij=dabs(epsij)
1190           aa_aq(i,j)=epsij*rrij*rrij
1191           bb_aq(i,j)=-sigeps*epsij*rrij
1192           aa_aq(j,i)=aa_aq(i,j)
1193           bb_aq(j,i)=bb_aq(i,j)
1194           epsijlip=epslip(i,j)
1195           sigeps=dsign(1.0D0,epsijlip)
1196           epsijlip=dabs(epsijlip)
1197           aa_lip(i,j)=epsijlip*rrij*rrij
1198           bb_lip(i,j)=-sigeps*epsijlip*rrij
1199           aa_lip(j,i)=aa_lip(i,j)
1200           bb_lip(j,i)=bb_lip(i,j)
1201           if (ipot.gt.2) then
1202             sigt1sq=sigma0(i)**2
1203             sigt2sq=sigma0(j)**2
1204             sigii1=sigii(i)
1205             sigii2=sigii(j)
1206             ratsig1=sigt2sq/sigt1sq
1207             ratsig2=1.0D0/ratsig1
1208             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1209             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1210             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1211           else
1212             rsum_max=sigma(i,j)
1213           endif
1214 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1215             sigmaii(i,j)=rsum_max
1216             sigmaii(j,i)=rsum_max
1217 c         else
1218 c           sigmaii(i,j)=r0(i,j)
1219 c           sigmaii(j,i)=r0(i,j)
1220 c         endif
1221 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1222           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1223             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1224             augm(i,j)=epsij*r_augm**(2*expon)
1225 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1226             augm(j,i)=augm(i,j)
1227           else
1228             augm(i,j)=0.0D0
1229             augm(j,i)=0.0D0
1230           endif
1231           if (lprint) then
1232             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1233      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1234      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1235           endif
1236         enddo
1237       enddo
1238 #ifdef OLDSCP
1239 C
1240 C Define the SC-p interaction constants (hard-coded; old style)
1241 C
1242       do i=1,ntyp
1243 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1244 C helix formation)
1245 c       aad(i,1)=0.3D0*4.0D0**12
1246 C Following line for constants currently implemented
1247 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1248         aad(i,1)=1.5D0*4.0D0**12
1249 c       aad(i,1)=0.17D0*5.6D0**12
1250         aad(i,2)=aad(i,1)
1251 C "Soft" SC-p repulsion
1252         bad(i,1)=0.0D0
1253 C Following line for constants currently implemented
1254 c       aad(i,1)=0.3D0*4.0D0**6
1255 C "Hard" SC-p repulsion
1256         bad(i,1)=3.0D0*4.0D0**6
1257 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1258         bad(i,2)=bad(i,1)
1259 c       aad(i,1)=0.0D0
1260 c       aad(i,2)=0.0D0
1261 c       bad(i,1)=1228.8D0
1262 c       bad(i,2)=1228.8D0
1263       enddo
1264 #else
1265 C
1266 C 8/9/01 Read the SC-p interaction constants from file
1267 C
1268       do i=1,ntyp
1269         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1270       enddo
1271       do i=1,ntyp
1272         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1273         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1274         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1275         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1276       enddo
1277 c      lprint=.true.
1278       if (lprint) then
1279         write (iout,*) "Parameters of SC-p interactions:"
1280         do i=1,ntyp
1281           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1282      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1283         enddo
1284       endif
1285 c      lprint=.false.
1286 #endif
1287 C
1288 C Define the constants of the disulfide bridge
1289 C
1290       ebr=-5.50D0
1291 c
1292 c Old arbitrary potential - commented out.
1293 c
1294 c      dbr= 4.20D0
1295 c      fbr= 3.30D0
1296 c
1297 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1298 c energy surface of diethyl disulfide.
1299 c A. Liwo and U. Kozlowska, 11/24/03
1300 c
1301       D0CM = 3.78d0
1302       AKCM = 15.1d0
1303       AKTH = 11.0d0
1304       AKCT = 12.0d0
1305       V1SS =-1.08d0
1306       V2SS = 7.61d0
1307       V3SS = 13.7d0
1308 c      akcm=0.0d0
1309 c      akth=0.0d0
1310 c      akct=0.0d0
1311 c      v1ss=0.0d0
1312 c      v2ss=0.0d0
1313 c      v3ss=0.0d0
1314
1315       if(me.eq.king) then
1316       write (iout,'(/a)') "Disulfide bridge parameters:"
1317       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1318       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1319       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1320       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1321      &  ' v3ss:',v3ss
1322       endif
1323       return
1324   111 write (iout,*) "Error reading bending energy parameters."
1325       goto 999
1326   112 write (iout,*) "Error reading rotamer energy parameters."
1327       goto 999
1328   113 write (iout,*) "Error reading torsional energy parameters."
1329       goto 999
1330   114 write (iout,*) "Error reading double torsional energy parameters."
1331       goto 999
1332   115 write (iout,*)
1333      &  "Error reading cumulant (multibody energy) parameters."
1334       goto 999
1335   116 write (iout,*) "Error reading electrostatic energy parameters."
1336       goto 999
1337  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1338       goto 999
1339   117 write (iout,*) "Error reading side chain interaction parameters."
1340       goto 999
1341   118 write (iout,*) "Error reading SCp interaction parameters."
1342       goto 999
1343   119 write (iout,*) "Error reading SCCOR parameters"
1344   999 continue
1345 #ifdef MPI
1346       call MPI_Finalize(Ierror)
1347 #endif
1348       stop
1349       return
1350       end
1351
1352
1353       subroutine getenv_loc(var, val)
1354       character(*) var, val
1355
1356 #ifdef WINIFL
1357       character(2000) line
1358       external ilen
1359
1360       open (196,file='env',status='old',readonly,shared)
1361       iread=0
1362 c      write(*,*)'looking for ',var
1363 10    read(196,*,err=11,end=11)line
1364       iread=index(line,var)
1365 c      write(*,*)iread,' ',var,' ',line
1366       if (iread.eq.0) go to 10
1367 c      write(*,*)'---> ',line
1368 11    continue
1369       if(iread.eq.0) then
1370 c       write(*,*)'CHUJ'
1371        val=''
1372       else
1373        iread=iread+ilen(var)+1
1374        read (line(iread:),*,err=12,end=12) val
1375 c       write(*,*)'OK: ',var,' = ',val
1376       endif
1377       close(196)
1378       return
1379 12    val=''
1380       close(196)
1381 #elif (defined CRAY)
1382       integer lennam,lenval,ierror
1383 c
1384 c        getenv using a POSIX call, useful on the T3D
1385 c        Sept 1996, comment out error check on advice of H. Pritchard
1386 c
1387       lennam = len(var)
1388       if(lennam.le.0) stop '--error calling getenv--'
1389       call pxfgetenv(var,lennam,val,lenval,ierror)
1390 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1391 #else
1392       call getenv(var,val)
1393 #endif
1394
1395       return
1396       end