source/wham/src-M-SAXS/parmread.F.safe

   1       subroutine parmread(iparm,*)
   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       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8       include 'DIMENSIONS.ZSCOPT'
   9       include 'DIMENSIONS.FREE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.CHAIN'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.GEO'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.TORSION'
  16       include 'COMMON.FFIELD'
  17       include 'COMMON.NAMES'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.WEIGHTS'
  20       include 'COMMON.ENEPS'
  21       include 'COMMON.SCCOR'
  22       include 'COMMON.SCROT'
  23       include 'COMMON.FREE'
  24       include 'COMMON.SHIELD'
  25       include 'COMMON.CONTROL'
  26       character*1 t1,t2,t3
  27       character*1 onelett(4) /"G","A","P","D"/
  28       character*1 toronelet(-2:2) /"p","a","G","A","P"/
  29       logical lprint
  30       dimension blower(3,3,maxlob)
  31       character*800 controlcard
  32       character*256 bondname_t,thetname_t,rotname_t,torname_t,
  33      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
  34      & sccorname_t
  35       integer ilen
  36       external ilen
  37       character*16 key
  38       integer iparm
  39       double precision ip,mp
  40       character*6 res1
  41       character*3 lancuch,ucase
  42 C      write (iout,*) "KURWA"
  43 C
  44 C Body
  45 C
  46       call getenv("PRINT_PARM",lancuch)
  47       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  48       write (iout,*) "lprint ",lprint
  49 C Set LPRINT=.TRUE. for debugging
  50       dwa16=2.0d0**(1.0d0/6.0d0)
  51       itypro=20
  52 C Assign virtual-bond length
  53       vbl=3.8D0
  54       vblinv=1.0D0/vbl
  55       vblinv2=vblinv*vblinv
  56       call card_concat(controlcard,.true.)
  57       wname(4)="WCORRH"
  58       do i=1,n_ene
  59         key = wname(i)(:ilen(wname(i)))
  60         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
  61       enddo
  62
  63       write (iout,*) "iparm",iparm," myparm",myparm
  64 c If reading not own parameters, skip assignment
  65       call reada(controlcard,"D0CM",d0cm,3.78d0)
  66       call reada(controlcard,"AKCM",akcm,15.1d0)
  67       call reada(controlcard,"AKTH",akth,11.0d0)
  68       call reada(controlcard,"AKCT",akct,12.0d0)
  69       call reada(controlcard,"V1SS",v1ss,-1.08d0)
  70       call reada(controlcard,"V2SS",v2ss,7.61d0)
  71       call reada(controlcard,"V3SS",v3ss,13.7d0)
  72       call reada(controlcard,"EBR",ebr,-5.50D0)
  73       call reada(controlcard,"DTRISS",dtriss,1.0D0)
  74       call reada(controlcard,"ATRISS",atriss,0.3D0)
  75       call reada(controlcard,"BTRISS",btriss,0.02D0)
  76       call reada(controlcard,"CTRISS",ctriss,1.0D0)
  77       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
  78       write(iout,*) "ATRISS",atriss
  79       write(iout,*) "BTRISS",btriss
  80       write(iout,*) "CTRISS",ctriss
  81       write(iout,*) "DTRISS",dtriss
  82
  83 C      do i=1,maxres
  84 C        dyn_ss_mask(i)=.false.
  85 C      enddo
  86 C      ebr=-12.0D0
  87 c
  88 c Old arbitrary potential - commented out.
  89 c
  90 c      dbr= 4.20D0
  91 c      fbr= 3.30D0
  92 c
  93 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
  94 c energy surface of diethyl disulfide.
  95 c A. Liwo and U. Kozlowska, 11/24/03
  96 c
  97       D0CM = 3.78d0
  98       AKCM = 15.1d0
  99       AKTH = 11.0d0
 100       AKCT = 12.0d0
 101       V1SS =-1.08d0
 102       V2SS = 7.61d0
 103       V3SS = 13.7d0
 104
 105       do i=1,maxres-1
 106         do j=i+1,maxres
 107           dyn_ssbond_ij(i,j)=1.0d300
 108         enddo
 109       enddo
 110       call reada(controlcard,"HT",Ht,0.0D0)
 111 C      if (dyn_ss) then
 112 C        ss_depth=ebr/wsc-0.25*eps(1,1)
 113 C        write(iout,*) HT,wsc,eps(1,1),'KURWA'
 114 C        Ht=Ht/wsc-0.25*eps(1,1)
 115
 116 C        akcm=akcm*whpb/wsc
 117 C        akth=akth*whpb/wsc
 118 C        akct=akct*whpb/wsc
 119 C        v1ss=v1ss*whpb/wsc
 120 C        v2ss=v2ss*whpb/wsc
 121 C        v3ss=v3ss*whpb/wsc
 122 C      else
 123 C        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
 124 C      endif
 125
 126       if (iparm.eq.myparm .or. .not.separate_parset) then
 127
 128 c
 129 c Setup weights for UNRES
 130 c
 131       wsc=ww(1)
 132       wscp=ww(2)
 133       welec=ww(3)
 134       wcorr=ww(4)
 135       wcorr5=ww(5)
 136       wcorr6=ww(6)
 137       wel_loc=ww(7)
 138       wturn3=ww(8)
 139       wturn4=ww(9)
 140       wturn6=ww(10)
 141       wang=ww(11)
 142       wscloc=ww(12)
 143       wtor=ww(13)
 144       wtor_d=ww(14)
 145       wvdwpp=ww(16)
 146       wbond=ww(18)
 147       wsccor=ww(19)
 148       whpb=ww(15)
 149       wstrain=ww(15)
 150       wliptran=ww(22)
 151       wshield=ww(25)
 152       endif
 153
 154       call card_concat(controlcard,.false.)
 155
 156 c Return if not own parameters
 157
 158       if (iparm.ne.myparm .and. separate_parset) return
 159
 160       call reads(controlcard,"BONDPAR",bondname_t,bondname)
 161       open (ibond,file=bondname_t,status='old')
 162       rewind(ibond)
 163       call reads(controlcard,"THETPAR",thetname_t,thetname)
 164       open (ithep,file=thetname_t,status='old')
 165       rewind(ithep)
 166       call reads(controlcard,"ROTPAR",rotname_t,rotname)
 167       open (irotam,file=rotname_t,status='old')
 168       rewind(irotam)
 169       call reads(controlcard,"TORPAR",torname_t,torname)
 170       open (itorp,file=torname_t,status='old')
 171       rewind(itorp)
 172       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 173       write (iout,*) "tor_mode",tor_mode
 174       call flush(iout)
 175       if (tor_mode.eq.0)
 176      & open (itordp,file=tordname_t,status='old')
 177       rewind(itordp)
 178       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
 179       open (isccor,file=sccorname_t,status='old')
 180       rewind(isccor)
 181       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 182       open (ifourier,file=fouriername_t,status='old')
 183       rewind(ifourier)
 184       call reads(controlcard,"ELEPAR",elename_t,elename)
 185       open (ielep,file=elename_t,status='old')
 186       rewind(ielep)
 187       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 188       open (isidep,file=sidename_t,status='old')
 189       rewind(isidep)
 190       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 191       open (iscpp,file=scpname_t,status='old')
 192       rewind(iscpp)
 193       write (iout,*) "Parameter set:",iparm
 194       write (iout,*) "Energy-term weights:"
 195       do i=1,n_ene
 196         write (iout,'(a16,f10.5)') wname(i),ww(i)
 197       enddo
 198       write (iout,*) "Sidechain potential file        : ",
 199      &  sidename_t(:ilen(sidename_t))
 200 #ifndef OLDSCP
 201       write (iout,*) "SCp potential file              : ",
 202      &  scpname_t(:ilen(scpname_t))
 203 #endif
 204       write (iout,*) "Electrostatic potential file    : ",
 205      &  elename_t(:ilen(elename_t))
 206       write (iout,*) "Cumulant coefficient file       : ",
 207      &  fouriername_t(:ilen(fouriername_t))
 208       write (iout,*) "Torsional parameter file        : ",
 209      &  torname_t(:ilen(torname_t))
 210       write (iout,*) "Double torsional parameter file : ",
 211      &  tordname_t(:ilen(tordname_t))
 212       write (iout,*) "Backbone-rotamer parameter file : ",
 213      &  sccorname(:ilen(sccorname))
 214       write (iout,*) "Bond & inertia constant file    : ",
 215      &  bondname_t(:ilen(bondname_t))
 216       write (iout,*) "Bending parameter file          : ",
 217      &  thetname_t(:ilen(thetname_t))
 218       write (iout,*) "Rotamer parameter file          : ",
 219      &  rotname_t(:ilen(rotname_t))
 220
 221 c
 222 c Read the virtual-bond parameters, masses, and moments of inertia
 223 c and Stokes' radii of the peptide group and side chains
 224 c
 225 #ifdef CRYST_BOND
 226       read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp
 227       do i=1,ntyp
 228         nbondterm(i)=1
 229         read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i)
 230         dsc(i) = vbldsc0(1,i)
 231         if (i.eq.10) then
 232           dsc_inv(i)=0.0D0
 233         else
 234           dsc_inv(i)=1.0D0/dsc(i)
 235         endif
 236       enddo
 237 #else
 238       read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk
 239       do i=1,ntyp
 240         read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
 241      &   aksc(j,i),abond0(j,i),j=1,nbondterm(i))
 242         dsc(i) = vbldsc0(1,i)
 243         if (i.eq.10) then
 244           dsc_inv(i)=0.0D0
 245         else
 246           dsc_inv(i)=1.0D0/dsc(i)
 247         endif
 248       enddo
 249 #endif
 250       if (lprint) then
 251         write(iout,'(/a/)')"Force constants virtual bonds:"
 252         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
 253      &   'inertia','Pstok'
 254         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 255         do i=1,ntyp
 256           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
 257      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
 258           do j=2,nbondterm(i)
 259             write (iout,'(13x,3f10.5)')
 260      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
 261           enddo
 262         enddo
 263       endif
 264        write (iout,*) "iliptranpar",iliptranpar
 265        write (iout,*) "liptranname ",liptranname
 266        read(iliptranpar,*,end=1161,err=1161) pepliptran
 267        write (iout,*) "pepliptran",pepliptran
 268        do i=1,ntyp
 269        read(iliptranpar,*,end=1161,err=1161) liptranene(i)
 270        write (iout,*) i,liptranene(i)
 271        enddo
 272        rewind iliptranpar
 273 #ifdef CRYST_THETA
 274 C
 275 C Read the parameters of the probability distribution/energy expression
 276 C of the virtual-bond valence angles theta
 277 C
 278       do i=1,ntyp
 279         read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
 280      &    (bthet(j,i,1,1),j=1,2)
 281         read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
 282         read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
 283         read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
 284         sigc0(i)=sigc0(i)**2
 285       enddo
 286       do i=1,ntyp
 287       athet(1,i,1,-1)=athet(1,i,1,1)
 288       athet(2,i,1,-1)=athet(2,i,1,1)
 289       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 290       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 291       athet(1,i,-1,1)=-athet(1,i,1,1)
 292       athet(2,i,-1,1)=-athet(2,i,1,1)
 293       bthet(1,i,-1,1)=bthet(1,i,1,1)
 294       bthet(2,i,-1,1)=bthet(2,i,1,1)
 295       enddo
 296       do i=-ntyp,-1
 297       a0thet(i)=a0thet(-i)
 298       athet(1,i,-1,-1)=athet(1,-i,1,1)
 299       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 300       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 301       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 302       athet(1,i,-1,1)=athet(1,-i,1,1)
 303       athet(2,i,-1,1)=-athet(2,-i,1,1)
 304       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 305       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 306       athet(1,i,1,-1)=-athet(1,-i,1,1)
 307       athet(2,i,1,-1)=athet(2,-i,1,1)
 308       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 309       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 310       theta0(i)=theta0(-i)
 311       sig0(i)=sig0(-i)
 312       sigc0(i)=sigc0(-i)
 313        do j=0,3
 314         polthet(j,i)=polthet(j,-i)
 315        enddo
 316        do j=1,3
 317          gthet(j,i)=gthet(j,-i)
 318        enddo
 319       enddo
 320       close (ithep)
 321       if (lprint) then
 322 c       write (iout,'(a)')
 323 c    &    'Parameters of the virtual-bond valence angles:'
 324 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 325 c    & '    ATHETA0   ','         A1   ','        A2    ',
 326 c    & '        B1    ','         B2   '
 327 c       do i=1,ntyp
 328 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 329 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 330 c       enddo
 331 c       write (iout,'(/a/9x,5a/79(1h-))')
 332 c    & 'Parameters of the expression for sigma(theta_c):',
 333 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 334 c    & '     ALPH3    ','    SIGMA0C   '
 335 c       do i=1,ntyp
 336 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 337 c    &      (polthet(j,i),j=0,3),sigc0(i)
 338 c       enddo
 339 c       write (iout,'(/a/9x,5a/79(1h-))')
 340 c    & 'Parameters of the second gaussian:',
 341 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 342 c    & '        G2    ','         G3   '
 343 c       do i=1,ntyp
 344 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 345 c    &       sig0(i),(gthet(j,i),j=1,3)
 346 c       enddo
 347         write (iout,'(a)')
 348      &    'Parameters of the virtual-bond valence angles:'
 349         write (iout,'(/a/9x,5a/79(1h-))')
 350      & 'Coefficients of expansion',
 351      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 352      & '   b1*10^1    ','    b2*10^1   '
 353         do i=1,ntyp
 354           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 355      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
 356      &        (10*bthet(j,i,1,1),j=1,2)
 357         enddo
 358         write (iout,'(/a/9x,5a/79(1h-))')
 359      & 'Parameters of the expression for sigma(theta_c):',
 360      & ' alpha0       ','  alph1       ',' alph2        ',
 361      & ' alhp3        ','   sigma0c    '
 362         do i=1,ntyp
 363           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 364      &      (polthet(j,i),j=0,3),sigc0(i)
 365         enddo
 366         write (iout,'(/a/9x,5a/79(1h-))')
 367      & 'Parameters of the second gaussian:',
 368      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 369      & '        G2    ','   G3*10^1    '
 370         do i=1,ntyp
 371           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 372      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 373         enddo
 374       endif
 375 #else
 376       IF (tor_mode.eq.0) THEN
 377 C
 378 C Read the parameters of Utheta determined from ab initio surfaces
 379 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 380 C
 381       read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
 382      &  ntheterm3,nsingle,ndouble
 383       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 384       read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
 385       do i=-ntyp1,-1
 386         ithetyp(i)=-ithetyp(-i)
 387       enddo
 388       do iblock=1,2
 389       do i=-maxthetyp,maxthetyp
 390         do j=-maxthetyp,maxthetyp
 391           do k=-maxthetyp,maxthetyp
 392             aa0thet(i,j,k,iblock)=0.0d0
 393             do l=1,ntheterm
 394               aathet(l,i,j,k,iblock)=0.0d0
 395             enddo
 396             do l=1,ntheterm2
 397               do m=1,nsingle
 398                 bbthet(m,l,i,j,k,iblock)=0.0d0
 399                 ccthet(m,l,i,j,k,iblock)=0.0d0
 400                 ddthet(m,l,i,j,k,iblock)=0.0d0
 401                 eethet(m,l,i,j,k,iblock)=0.0d0
 402               enddo
 403             enddo
 404             do l=1,ntheterm3
 405               do m=1,ndouble
 406                 do mm=1,ndouble
 407                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 408                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 409                 enddo
 410               enddo
 411             enddo
 412           enddo
 413         enddo
 414       enddo
 415       enddo
 416 C      write (iout,*) "KURWA1"
 417       do iblock=1,2
 418       do i=0,nthetyp
 419         do j=-nthetyp,nthetyp
 420           do k=-nthetyp,nthetyp
 421             read (ithep,'(6a)',end=111,err=111) res1
 422             write(iout,*) res1,i,j,k
 423             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
 424             read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
 425      &        l=1,ntheterm)
 426             read (ithep,*,end=111,err=111)
 427      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 428      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 429      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 430      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 431      &        ,ll=1,ntheterm2)
 432             read (ithep,*,end=111,err=111)
 433      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 434      &      ffthet(lll,llll,ll,i,j,k,iblock),
 435      &         ggthet(llll,lll,ll,i,j,k,iblock)
 436      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 437      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 438           enddo
 439         enddo
 440       enddo
 441 C       write(iout,*) "KURWA1.1"
 442 C
 443 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 444 C coefficients of theta-and-gamma-dependent terms are zero.
 445 C
 446       do i=1,nthetyp
 447         do j=1,nthetyp
 448           do l=1,ntheterm
 449             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 450             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 451           enddo
 452           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 453           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 454         enddo
 455         do l=1,ntheterm
 456           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 457         enddo
 458         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 459       enddo
 460       enddo
 461 C       write(iout,*) "KURWA1.5"
 462 C Substitution for D aminoacids from symmetry.
 463       do iblock=1,2
 464       do i=-nthetyp,0
 465         do j=-nthetyp,nthetyp
 466           do k=-nthetyp,nthetyp
 467            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 468            do l=1,ntheterm
 469            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 470            enddo
 471            do ll=1,ntheterm2
 472             do lll=1,nsingle
 473             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 474             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 475             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 476             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 477             enddo
 478           enddo
 479           do ll=1,ntheterm3
 480            do lll=2,ndouble
 481             do llll=1,lll-1
 482             ffthet(llll,lll,ll,i,j,k,iblock)=
 483      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 484             ffthet(lll,llll,ll,i,j,k,iblock)=
 485      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 486             ggthet(llll,lll,ll,i,j,k,iblock)=
 487      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 488             ggthet(lll,llll,ll,i,j,k,iblock)=
 489      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 490             enddo !ll
 491            enddo  !lll
 492           enddo   !llll
 493          enddo    !k
 494         enddo     !j
 495        enddo      !i
 496       enddo       !iblock
 497
 498 C
 499 C Control printout of the coefficients of virtual-bond-angle potentials
 500 C
 501       if (lprint) then
 502         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 503         do iblock=1,2
 504         do i=1,nthetyp+1
 505           do j=1,nthetyp+1
 506             do k=1,nthetyp+1
 507               write (iout,'(//4a)')
 508      &         'Type ',onelett(i),onelett(j),onelett(k)
 509               write (iout,'(//a,10x,a)') " l","a[l]"
 510               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 511               write (iout,'(i2,1pe15.5)')
 512      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 513             do l=1,ntheterm2
 514               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 515      &          "b",l,"c",l,"d",l,"e",l
 516               do m=1,nsingle
 517                 write (iout,'(i2,4(1pe15.5))') m,
 518      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 519      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 520               enddo
 521             enddo
 522             do l=1,ntheterm3
 523               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 524      &          "f+",l,"f-",l,"g+",l,"g-",l
 525               do m=2,ndouble
 526                 do n=1,m-1
 527                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 528      &              ffthet(n,m,l,i,j,k,iblock),
 529      &              ffthet(m,n,l,i,j,k,iblock),
 530      &              ggthet(n,m,l,i,j,k,iblock),
 531      &              ggthet(m,n,l,i,j,k,iblock)
 532                 enddo
 533               enddo
 534             enddo
 535           enddo
 536         enddo
 537       enddo
 538       enddo
 539       call flush(iout)
 540       endif
 541
 542       ELSE
 543
 544 C here will be the apropriate recalibrating for D-aminoacid
 545       read (ithep,*,end=111,err=111) nthetyp
 546       do i=-nthetyp+1,nthetyp-1
 547         read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
 548         do j=0,nbend_kcc_Tb(i)
 549           read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
 550         enddo
 551       enddo
 552       if (lprint) then
 553         write (iout,'(a)')
 554      &    "Parameters of the valence-only potentials"
 555         do i=-nthetyp+1,nthetyp-1
 556         write (iout,'(2a)') "Type ",toronelet(i)
 557         do k=0,nbend_kcc_Tb(i)
 558           write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
 559         enddo
 560         enddo
 561       endif
 562
 563       ENDIF ! TOR_MODE
 564
 565       close(ithep)
 566 #endif
 567 C      write(iout,*) 'KURWA2'
 568 #ifdef CRYST_SC
 569 C
 570 C Read the parameters of the probability distribution/energy expression
 571 C of the side chains.
 572 C
 573       do i=1,ntyp
 574 cc      write (iout,*) "tu dochodze",i
 575         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 576         if (i.eq.10) then
 577           dsc_inv(i)=0.0D0
 578         else
 579           dsc_inv(i)=1.0D0/dsc(i)
 580         endif
 581         if (i.ne.10) then
 582         do j=1,nlob(i)
 583           do k=1,3
 584             do l=1,3
 585               blower(l,k,j)=0.0D0
 586             enddo
 587           enddo
 588         enddo
 589         bsc(1,i)=0.0D0
 590         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 591      &    ((blower(k,l,1),l=1,k),k=1,3)
 592         censc(1,1,-i)=censc(1,1,i)
 593         censc(2,1,-i)=censc(2,1,i)
 594         censc(3,1,-i)=-censc(3,1,i)
 595         do j=2,nlob(i)
 596           read (irotam,*,end=112,err=112) bsc(j,i)
 597           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 598      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 599         censc(1,j,-i)=censc(1,j,i)
 600         censc(2,j,-i)=censc(2,j,i)
 601         censc(3,j,-i)=-censc(3,j,i)
 602 C BSC is amplitude of Gaussian
 603         enddo
 604         do j=1,nlob(i)
 605           do k=1,3
 606             do l=1,k
 607               akl=0.0D0
 608               do m=1,3
 609                 akl=akl+blower(k,m,j)*blower(l,m,j)
 610               enddo
 611               gaussc(k,l,j,i)=akl
 612               gaussc(l,k,j,i)=akl
 613              if (((k.eq.3).and.(l.ne.3))
 614      &        .or.((l.eq.3).and.(k.ne.3))) then
 615                 gaussc(k,l,j,-i)=-akl
 616                 gaussc(l,k,j,-i)=-akl
 617               else
 618                 gaussc(k,l,j,-i)=akl
 619                 gaussc(l,k,j,-i)=akl
 620               endif
 621             enddo
 622           enddo
 623         enddo
 624         endif
 625       enddo
 626       close (irotam)
 627       if (lprint) then
 628         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 629         do i=1,ntyp
 630           nlobi=nlob(i)
 631           if (nlobi.gt.0) then
 632           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 633      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 634 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 635 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 636 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 637 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 638            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 639      &                             'log h',(bsc(j,i),j=1,nlobi)
 640            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 641      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 642 c          write (iout,'(a)')
 643 c         do j=1,nlobi
 644 c           ind=0
 645 c           do k=1,3
 646 c             do l=1,k
 647 c              ind=ind+1
 648 c              blower(k,l,j)=gaussc(ind,j,i)
 649 c             enddo
 650 c           enddo
 651 c         enddo
 652           do k=1,3
 653             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 654      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 655           enddo
 656           endif
 657         enddo
 658       endif
 659 #else
 660 C
 661 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 662 C added by Urszula Kozlowska 07/11/2007
 663 C
 664       do i=1,ntyp
 665         read (irotam,*,end=112,err=112)
 666        if (i.eq.10) then
 667          read (irotam,*,end=112,err=112)
 668        else
 669          do j=1,65
 670            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 671          enddo
 672        endif
 673       enddo
 674 #endif
 675       close(irotam)
 676 C
 677 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 678 C         interaction energy of the Gly, Ala, and Pro prototypes.
 679 C
 680       read (ifourier,*,end=115,err=115) nloctyp
 681 #ifdef NEWCORR
 682       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
 683       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
 684       itype2loc(ntyp1)=nloctyp
 685       iloctyp(nloctyp)=ntyp1
 686       do i=1,ntyp1
 687         itype2loc(-i)=-itype2loc(i)
 688       enddo
 689 #else
 690       iloctyp(0)=10
 691       iloctyp(1)=9
 692       iloctyp(2)=20
 693       iloctyp(3)=ntyp1
 694 #endif
 695       do i=1,nloctyp
 696         iloctyp(-i)=-iloctyp(i)
 697       enddo
 698 c      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
 699 c      write (iout,*) "nloctyp",nloctyp,
 700 c     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
 701 #ifdef NEWCORR
 702       do i=0,nloctyp-1
 703 c             write (iout,*) "NEWCORR",i
 704         read (ifourier,*,end=115,err=115)
 705         do ii=1,3
 706           do j=1,2
 707             read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
 708           enddo
 709         enddo
 710 c             write (iout,*) "NEWCORR BNEW1"
 711 c             write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
 712         do ii=1,3
 713           do j=1,2
 714             read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
 715           enddo
 716         enddo
 717 c             write (iout,*) "NEWCORR BNEW2"
 718 c             write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
 719         do kk=1,3
 720           read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
 721           read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
 722         enddo
 723 c             write (iout,*) "NEWCORR CCNEW"
 724 c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
 725         do kk=1,3
 726           read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
 727           read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
 728         enddo
 729 c             write (iout,*) "NEWCORR DDNEW"
 730 c             write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
 731         do ii=1,2
 732           do jj=1,2
 733             do kk=1,2
 734               read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
 735             enddo
 736           enddo
 737         enddo
 738 c             write (iout,*) "NEWCORR EENEW1"
 739 c             write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
 740         do ii=1,3
 741           read (ifourier,*,end=115,err=115) e0new(ii,i)
 742         enddo
 743 c             write (iout,*) (e0new(ii,i),ii=1,3)
 744       enddo
 745 c             write (iout,*) "NEWCORR EENEW"
 746       do i=0,nloctyp-1
 747         do ii=1,3
 748           ccnew(ii,1,i)=ccnew(ii,1,i)/2
 749           ccnew(ii,2,i)=ccnew(ii,2,i)/2
 750           ddnew(ii,1,i)=ddnew(ii,1,i)/2
 751           ddnew(ii,2,i)=ddnew(ii,2,i)/2
 752         enddo
 753       enddo
 754       do i=1,nloctyp-1
 755         do ii=1,3
 756           bnew1(ii,1,-i)= bnew1(ii,1,i)
 757           bnew1(ii,2,-i)=-bnew1(ii,2,i)
 758           bnew2(ii,1,-i)= bnew2(ii,1,i)
 759           bnew2(ii,2,-i)=-bnew2(ii,2,i)
 760         enddo
 761         do ii=1,3
 762 c          ccnew(ii,1,i)=ccnew(ii,1,i)/2
 763 c          ccnew(ii,2,i)=ccnew(ii,2,i)/2
 764 c          ddnew(ii,1,i)=ddnew(ii,1,i)/2
 765 c          ddnew(ii,2,i)=ddnew(ii,2,i)/2
 766           ccnew(ii,1,-i)=ccnew(ii,1,i)
 767           ccnew(ii,2,-i)=-ccnew(ii,2,i)
 768           ddnew(ii,1,-i)=ddnew(ii,1,i)
 769           ddnew(ii,2,-i)=-ddnew(ii,2,i)
 770         enddo
 771         e0new(1,-i)= e0new(1,i)
 772         e0new(2,-i)=-e0new(2,i)
 773         e0new(3,-i)=-e0new(3,i)
 774         do kk=1,2
 775           eenew(kk,1,1,-i)= eenew(kk,1,1,i)
 776           eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
 777           eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
 778           eenew(kk,2,2,-i)= eenew(kk,2,2,i)
 779         enddo
 780       enddo
 781       if (lprint) then
 782         write (iout,'(a)') "Coefficients of the multibody terms"
 783         do i=-nloctyp+1,nloctyp-1
 784           write (iout,*) "Type: ",onelet(iloctyp(i))
 785           write (iout,*) "Coefficients of the expansion of B1"
 786           do j=1,2
 787             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
 788           enddo
 789           write (iout,*) "Coefficients of the expansion of B2"
 790           do j=1,2
 791             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
 792           enddo
 793           write (iout,*) "Coefficients of the expansion of C"
 794           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
 795           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
 796           write (iout,*) "Coefficients of the expansion of D"
 797           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
 798           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
 799           write (iout,*) "Coefficients of the expansion of E"
 800           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
 801           do j=1,2
 802             do k=1,2
 803               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
 804             enddo
 805           enddo
 806         enddo
 807       endif
 808 #else
 809       if (lprint)
 810      &  write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
 811       do i=0,nloctyp-1
 812         read (ifourier,*,end=115,err=115)
 813         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
 814         if (lprint) then
 815         write (iout,*) 'Type ',onelet(iloctyp(i))
 816         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 817         endif
 818         if (i.gt.0) then
 819         b(2,-i)= b(2,i)
 820         b(3,-i)= b(3,i)
 821         b(4,-i)=-b(4,i)
 822         b(5,-i)=-b(5,i)
 823         endif
 824 c        B1(1,i)  = b(3)
 825 c        B1(2,i)  = b(5)
 826 c        B1(1,-i) = b(3)
 827 c        B1(2,-i) = -b(5)
 828 c        b1(1,i)=0.0d0
 829 c        b1(2,i)=0.0d0
 830 c        B1tilde(1,i) = b(3)
 831 c        B1tilde(2,i) =-b(5)
 832 c        B1tilde(1,-i) =-b(3)
 833 c        B1tilde(2,-i) =b(5)
 834 c        b1tilde(1,i)=0.0d0
 835 c        b1tilde(2,i)=0.0d0
 836 c        B2(1,i)  = b(2)
 837 c        B2(2,i)  = b(4)
 838 c        B2(1,-i)  =b(2)
 839 c        B2(2,-i)  =-b(4)
 840 cc        B1tilde(1,i) = b(3,i)
 841 cc        B1tilde(2,i) =-b(5,i)
 842 C        B1tilde(1,-i) =-b(3,i)
 843 C        B1tilde(2,-i) =b(5,i)
 844 cc        b1tilde(1,i)=0.0d0
 845 cc        b1tilde(2,i)=0.0d0
 846 cc        B2(1,i)  = b(2,i)
 847 cc        B2(2,i)  = b(4,i)
 848 C        B2(1,-i)  =b(2,i)
 849 C        B2(2,-i)  =-b(4,i)
 850
 851 c        b2(1,i)=0.0d0
 852 c        b2(2,i)=0.0d0
 853         CCold(1,1,i)= b(7,i)
 854         CCold(2,2,i)=-b(7,i)
 855         CCold(2,1,i)= b(9,i)
 856         CCold(1,2,i)= b(9,i)
 857         CCold(1,1,-i)= b(7,i)
 858         CCold(2,2,-i)=-b(7,i)
 859         CCold(2,1,-i)=-b(9,i)
 860         CCold(1,2,-i)=-b(9,i)
 861 c        CC(1,1,i)=0.0d0
 862 c        CC(2,2,i)=0.0d0
 863 c        CC(2,1,i)=0.0d0
 864 c        CC(1,2,i)=0.0d0
 865 c        Ctilde(1,1,i)= CCold(1,1,i)
 866 c        Ctilde(1,2,i)= CCold(1,2,i)
 867 c        Ctilde(2,1,i)=-CCold(2,1,i)
 868 c        Ctilde(2,2,i)=-CCold(2,2,i)
 869
 870 c        Ctilde(1,1,i)=0.0d0
 871 c        Ctilde(1,2,i)=0.0d0
 872 c        Ctilde(2,1,i)=0.0d0
 873 c        Ctilde(2,2,i)=0.0d0
 874         DDold(1,1,i)= b(6,i)
 875         DDold(2,2,i)=-b(6,i)
 876         DDold(2,1,i)= b(8,i)
 877         DDold(1,2,i)= b(8,i)
 878         DDold(1,1,-i)= b(6,i)
 879         DDold(2,2,-i)=-b(6,i)
 880         DDold(2,1,-i)=-b(8,i)
 881         DDold(1,2,-i)=-b(8,i)
 882 c        DD(1,1,i)=0.0d0
 883 c        DD(2,2,i)=0.0d0
 884 c        DD(2,1,i)=0.0d0
 885 c        DD(1,2,i)=0.0d0
 886 c        Dtilde(1,1,i)= DD(1,1,i)
 887 c        Dtilde(1,2,i)= DD(1,2,i)
 888 c        Dtilde(2,1,i)=-DD(2,1,i)
 889 c        Dtilde(2,2,i)=-DD(2,2,i)
 890
 891 c        Dtilde(1,1,i)=0.0d0
 892 c        Dtilde(1,2,i)=0.0d0
 893 c        Dtilde(2,1,i)=0.0d0
 894 c        Dtilde(2,2,i)=0.0d0
 895         EEold(1,1,i)= b(10,i)+b(11,i)
 896         EEold(2,2,i)=-b(10,i)+b(11,i)
 897         EEold(2,1,i)= b(12,i)-b(13,i)
 898         EEold(1,2,i)= b(12,i)+b(13,i)
 899         EEold(1,1,-i)= b(10,i)+b(11,i)
 900         EEold(2,2,-i)=-b(10,i)+b(11,i)
 901         EEold(2,1,-i)=-b(12,i)+b(13,i)
 902         EEold(1,2,-i)=-b(12,i)-b(13,i)
 903         write(iout,*) "TU DOCHODZE"
 904         print *,"JESTEM"
 905 c        ee(1,1,i)=1.0d0
 906 c        ee(2,2,i)=1.0d0
 907 c        ee(2,1,i)=0.0d0
 908 c        ee(1,2,i)=0.0d0
 909 c        ee(2,1,i)=ee(1,2,i)
 910       enddo
 911       if (lprint) then
 912       write (iout,*)
 913       write (iout,*)
 914      &"Coefficients of the cumulants (independent of valence angles)"
 915       do i=-nloctyp+1,nloctyp-1
 916         write (iout,*) 'Type ',onelet(iloctyp(i))
 917         write (iout,*) 'B1'
 918         write(iout,'(2f10.5)') B(3,i),B(5,i)
 919         write (iout,*) 'B2'
 920         write(iout,'(2f10.5)') B(2,i),B(4,i)
 921         write (iout,*) 'CC'
 922         do j=1,2
 923           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
 924         enddo
 925         write(iout,*) 'DD'
 926         do j=1,2
 927           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
 928         enddo
 929         write(iout,*) 'EE'
 930         do j=1,2
 931           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
 932         enddo
 933       enddo
 934       endif
 935 #endif
 936 C      write (iout,*) 'KURWAKURWA'
 937 #ifdef CRYST_TOR
 938 C
 939 C Read torsional parameters in old format
 940 C
 941       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 942       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 943       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 944       do i=1,ntortyp
 945         do j=1,ntortyp
 946           read (itorp,'(a)',end=113,err=113)
 947           do k=1,nterm_old
 948             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 949           enddo
 950         enddo
 951       enddo
 952       close (itorp)
 953       if (lprint) then
 954         write (iout,'(/a/)') 'Torsional constants:'
 955         do i=1,ntortyp
 956           do j=1,ntortyp
 957             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 958             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 959           enddo
 960         enddo
 961       endif
 962 #else
 963 C
 964 C Read torsional parameters
 965 C
 966       IF (TOR_MODE.eq.0) THEN
 967
 968       read (itorp,*,end=113,err=113) ntortyp
 969       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 970       do i=1,ntyp1
 971         itype2loc(i)=itortyp(i)
 972       enddo
 973       write (iout,*) 'ntortyp',ntortyp
 974       do i=1,ntyp1
 975         itype2loc(-i)=-itype2loc(i)
 976       enddo
 977       itortyp(ntyp1)=ntortyp
 978       do iblock=1,2
 979       do i=-ntyp,-1
 980        itortyp(i)=-itortyp(-i)
 981       enddo
 982 c      write (iout,*) 'ntortyp',ntortyp
 983       do i=0,ntortyp-1
 984         do j=-ntortyp+1,ntortyp-1
 985           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 986      &          nlor(i,j,iblock)
 987           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 988           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 989           v0ij=0.0d0
 990           si=-1.0d0
 991           do k=1,nterm(i,j,iblock)
 992             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 993      &      v2(k,i,j,iblock)
 994             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 995             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 996             v0ij=v0ij+si*v1(k,i,j,iblock)
 997             si=-si
 998          enddo
 999           do k=1,nlor(i,j,iblock)
1000             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1001      &        vlor2(k,i,j),vlor3(k,i,j)
1002             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1003           enddo
1004           v0(i,j,iblock)=v0ij
1005           v0(-i,-j,iblock)=v0ij
1006         enddo
1007       enddo
1008       enddo
1009       close (itorp)
1010       if (lprint) then
1011         write (iout,'(/a/)') 'Torsional constants:'
1012         do i=1,ntortyp
1013           do j=1,ntortyp
1014             write (iout,*) 'ityp',i,' jtyp',j
1015             write (iout,*) 'Fourier constants'
1016             do k=1,nterm(i,j,iblock)
1017               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1018      &        v2(k,i,j,iblock)
1019             enddo
1020             write (iout,*) 'Lorenz constants'
1021             do k=1,nlor(i,j,iblock)
1022               write (iout,'(3(1pe15.5))')
1023      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1024             enddo
1025           enddo
1026         enddo
1027       endif
1028 C
1029 C 6/23/01 Read parameters for double torsionals
1030 C
1031       do iblock=1,2
1032       do i=0,ntortyp-1
1033         do j=-ntortyp+1,ntortyp-1
1034           do k=-ntortyp+1,ntortyp-1
1035             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1036 c              write (iout,*) "OK onelett",
1037 c     &         i,j,k,t1,t2,t3
1038
1039             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1040      &        .or. t3.ne.toronelet(k)) then
1041               write (iout,*) "Error in double torsional parameter file",
1042      &         i,j,k,t1,t2,t3
1043 #ifdef MPI
1044               call MPI_Finalize(Ierror)
1045 #endif
1046                stop "Error in double torsional parameter file"
1047             endif
1048           read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1049      &         ntermd_2(i,j,k,iblock)
1050             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1051             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1052             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1053      &         ntermd_1(i,j,k,iblock))
1054             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1055      &         ntermd_1(i,j,k,iblock))
1056             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1057      &         ntermd_1(i,j,k,iblock))
1058             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1059      &         ntermd_1(i,j,k,iblock))
1060 C Martix of D parameters for one dimesional foureir series
1061             do l=1,ntermd_1(i,j,k,iblock)
1062              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1063              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1064              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1065              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1066 c            write(iout,*) "whcodze" ,
1067 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1068             enddo
1069             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1070      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1071      &         v2s(m,l,i,j,k,iblock),
1072      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1073 C Martix of D parameters for two dimesional fourier series
1074             do l=1,ntermd_2(i,j,k,iblock)
1075              do m=1,l-1
1076              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1077              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1078              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1079              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1080              enddo!m
1081             enddo!l
1082           enddo!k
1083         enddo!j
1084       enddo!i
1085       enddo!iblock
1086       if (lprint) then
1087       write (iout,*)
1088       write (iout,*) 'Constants for double torsionals'
1089       do iblock=1,2
1090       do i=0,ntortyp-1
1091         do j=-ntortyp+1,ntortyp-1
1092           do k=-ntortyp+1,ntortyp-1
1093             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1094      &        ' nsingle',ntermd_1(i,j,k,iblock),
1095      &        ' ndouble',ntermd_2(i,j,k,iblock)
1096             write (iout,*)
1097             write (iout,*) 'Single angles:'
1098             do l=1,ntermd_1(i,j,k,iblock)
1099               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1100      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1101      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1102      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1103             enddo
1104             write (iout,*)
1105             write (iout,*) 'Pairs of angles:'
1106             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1107             do l=1,ntermd_2(i,j,k,iblock)
1108               write (iout,'(i5,20f10.5)')
1109      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1110             enddo
1111             write (iout,*)
1112            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1113             do l=1,ntermd_2(i,j,k,iblock)
1114               write (iout,'(i5,20f10.5)')
1115      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1116      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1117             enddo
1118             write (iout,*)
1119           enddo
1120         enddo
1121       enddo
1122       enddo
1123       endif
1124
1125       ELSE IF (TOR_MODE.eq.1) THEN
1126
1127 C read valence-torsional parameters
1128       read (itorp,*,end=113,err=113) ntortyp
1129       nkcctyp=ntortyp
1130       write (iout,*) "Valence-torsional parameters read in ntortyp",
1131      &   ntortyp
1132       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1133       write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1134       do i=-ntyp,-1
1135         itortyp(i)=-itortyp(-i)
1136       enddo
1137       do i=-ntortyp+1,ntortyp-1
1138         do j=-ntortyp+1,ntortyp-1
1139 C first we read the cos and sin gamma parameters
1140           read (itorp,'(13x,a)',end=113,err=113) string
1141           write (iout,*) i,j,string
1142           read (itorp,*,end=113,err=113)
1143      &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1144 C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1145           do k=1,nterm_kcc(j,i)
1146             do l=1,nterm_kcc_Tb(j,i)
1147               do ll=1,nterm_kcc_Tb(j,i)
1148               read (itorp,*,end=113,err=113) ii,jj,kk,
1149      &          v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1150               enddo
1151             enddo
1152           enddo
1153         enddo
1154       enddo
1155
1156       ELSE
1157 c AL 4/8/16: Calculate coefficients from one-body parameters
1158       ntortyp=nloctyp
1159       do i=-ntyp1,ntyp1
1160         itortyp(i)=itype2loc(i)
1161       enddo
1162       write (iout,*)
1163      &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1164      & ,ntortyp
1165       do i=-ntortyp+1,ntortyp-1
1166         do j=-ntortyp+1,ntortyp-1
1167           nterm_kcc(j,i)=2
1168           nterm_kcc_Tb(j,i)=3
1169           do k=1,nterm_kcc_Tb(j,i)
1170             do l=1,nterm_kcc_Tb(j,i)
1171               v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j)
1172      &                         +bnew1(k,2,i)*bnew2(l,2,j)
1173               v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j)
1174      &                         +bnew1(k,2,i)*bnew2(l,1,j)
1175             enddo
1176           enddo
1177           do k=1,nterm_kcc_Tb(j,i)
1178             do l=1,nterm_kcc_Tb(j,i)
1179               v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j)
1180      &                         -ccnew(k,2,i)*ddnew(l,2,j))
1181               v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j)
1182      &                         +ccnew(k,1,i)*ddnew(l,2,j))
1183             enddo
1184           enddo
1185 cf(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1186         enddo
1187       enddo
1188
1189       ENDIF ! TOR_MODE
1190
1191       if (tor_mode.gt.0 .and. lprint) then
1192 c Print valence-torsional parameters
1193         write (iout,'(a)')
1194      &    "Parameters of the valence-torsional potentials"
1195         do i=-ntortyp+1,ntortyp-1
1196         do j=-ntortyp+1,ntortyp-1
1197         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1198         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1199         do k=1,nterm_kcc(j,i)
1200           do l=1,nterm_kcc_Tb(j,i)
1201             do ll=1,nterm_kcc_Tb(j,i)
1202                write (iout,'(3i5,2f15.4)')
1203      &            k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1204             enddo
1205           enddo
1206         enddo
1207         enddo
1208         enddo
1209       endif
1210
1211 #endif
1212 C Read of Side-chain backbone correlation parameters
1213 C Modified 11 May 2012 by Adasko
1214 CCC
1215 C
1216       read (isccor,*,end=119,err=119) nsccortyp
1217       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1218       do i=-ntyp,-1
1219         isccortyp(i)=-isccortyp(-i)
1220       enddo
1221       iscprol=isccortyp(20)
1222 c      write (iout,*) 'ntortyp',ntortyp
1223       maxinter=3
1224 cc maxinter is maximum interaction sites
1225       do l=1,maxinter
1226       do i=1,nsccortyp
1227         do j=1,nsccortyp
1228           read (isccor,*,end=119,err=119)
1229      &nterm_sccor(i,j),nlor_sccor(i,j)
1230 c          write (iout,*) nterm_sccor(i,j)
1231           v0ijsccor=0.0d0
1232           v0ijsccor1=0.0d0
1233           v0ijsccor2=0.0d0
1234           v0ijsccor3=0.0d0
1235           si=-1.0d0
1236           nterm_sccor(-i,j)=nterm_sccor(i,j)
1237           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1238           nterm_sccor(i,-j)=nterm_sccor(i,j)
1239 c          write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1240 c     &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
1241           do k=1,nterm_sccor(i,j)
1242             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1243      &    ,v2sccor(k,l,i,j)
1244             if (j.eq.iscprol) then
1245              if (i.eq.isccortyp(10)) then
1246              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1247              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1248              else
1249              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1250      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1251              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1252      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1253              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1254              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1255              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1256              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1257              endif
1258             else
1259              if (i.eq.isccortyp(10)) then
1260              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1261              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1262              else
1263                if (j.eq.isccortyp(10)) then
1264              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1265              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1266                else
1267              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1268              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1269              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1270              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1271              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1272              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1273                 endif
1274                endif
1275             endif
1276             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1277             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1278             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1279             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1280             si=-si
1281            enddo
1282           do k=1,nlor_sccor(i,j)
1283             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1284      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1285             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1286      &(1+vlor3sccor(k,i,j)**2)
1287           enddo
1288           v0sccor(l,i,j)=v0ijsccor
1289           v0sccor(l,-i,j)=v0ijsccor1
1290           v0sccor(l,i,-j)=v0ijsccor2
1291           v0sccor(l,-i,-j)=v0ijsccor3
1292           enddo
1293         enddo
1294       enddo
1295       close (isccor)
1296       if (lprint) then
1297         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1298         do l=1,maxinter
1299         do i=1,nsccortyp
1300           do j=1,nsccortyp
1301             write (iout,*) 'ityp',i,' jtyp',j
1302             write (iout,*) 'Fourier constants'
1303             do k=1,nterm_sccor(i,j)
1304               write (iout,'(2(1pe15.5))')
1305      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1306             enddo
1307             write (iout,*) 'Lorenz constants'
1308             do k=1,nlor_sccor(i,j)
1309               write (iout,'(3(1pe15.5))')
1310      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1311             enddo
1312           enddo
1313         enddo
1314         enddo
1315       endif
1316 C
1317 C Read electrostatic-interaction parameters
1318 C
1319       if (lprint) then
1320         write (iout,'(/a)') 'Electrostatic interaction constants:'
1321         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1322      &            'IT','JT','APP','BPP','AEL6','AEL3'
1323       endif
1324       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1325       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1326       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1327       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1328       close (ielep)
1329       do i=1,2
1330         do j=1,2
1331         rri=rpp(i,j)**6
1332         app (i,j)=epp(i,j)*rri*rri
1333         bpp (i,j)=-2.0D0*epp(i,j)*rri
1334         ael6(i,j)=elpp6(i,j)*4.2D0**6
1335         ael3(i,j)=elpp3(i,j)*4.2D0**3
1336         lprint=.true.
1337         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1338      &                    ael6(i,j),ael3(i,j)
1339          lprint=.false.
1340         enddo
1341       enddo
1342 C
1343 C Read side-chain interaction parameters.
1344 C
1345       read (isidep,*,end=117,err=117) ipot,expon
1346       if (ipot.lt.1 .or. ipot.gt.5) then
1347         write (iout,'(2a)') 'Error while reading SC interaction',
1348      &               'potential file - unknown potential type.'
1349         stop
1350       endif
1351       expon2=expon/2
1352       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1353      & ', exponents are ',expon,2*expon
1354       goto (10,20,30,30,40) ipot
1355 C----------------------- LJ potential ---------------------------------
1356    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1357      &   (sigma0(i),i=1,ntyp)
1358       if (lprint) then
1359         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1360         write (iout,'(a/)') 'The epsilon array:'
1361         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1362         write (iout,'(/a)') 'One-body parameters:'
1363         write (iout,'(a,4x,a)') 'residue','sigma'
1364         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1365       endif
1366       goto 50
1367 C----------------------- LJK potential --------------------------------
1368    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1369      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1370       if (lprint) then
1371         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1372         write (iout,'(a/)') 'The epsilon array:'
1373         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1374         write (iout,'(/a)') 'One-body parameters:'
1375         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1376         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1377      &        i=1,ntyp)
1378       endif
1379       goto 50
1380 C---------------------- GB or BP potential -----------------------------
1381    30 do i=1,ntyp
1382        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1383       enddo
1384       read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1385       read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1386       read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1387       read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1388       do i=1,ntyp
1389        read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1390 C       write(iout,*) "WARNING!!",i,ntyp
1391        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1392 C       do j=1,ntyp
1393 C       epslip(i,j)=epslip(i,j)+0.05d0
1394 C       enddo
1395       enddo
1396 C For the GB potential convert sigma'**2 into chi'
1397       if (ipot.eq.4) then
1398         do i=1,ntyp
1399           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1400         enddo
1401       endif
1402       if (lprint) then
1403         write (iout,'(/a/)') 'Parameters of the BP potential:'
1404         write (iout,'(a/)') 'The epsilon array:'
1405         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1406         write (iout,'(/a)') 'One-body parameters:'
1407         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1408      &       '    chip  ','    alph  '
1409         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1410      &                     chip(i),alp(i),i=1,ntyp)
1411       endif
1412       goto 50
1413 C--------------------- GBV potential -----------------------------------
1414    40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1415      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1416      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1417       if (lprint) then
1418         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1419         write (iout,'(a/)') 'The epsilon array:'
1420         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1421         write (iout,'(/a)') 'One-body parameters:'
1422         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1423      &      's||/s_|_^2','    chip  ','    alph  '
1424         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1425      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1426       endif
1427    50 continue
1428       close (isidep)
1429 C-----------------------------------------------------------------------
1430 C Calculate the "working" parameters of SC interactions.
1431       do i=2,ntyp
1432         do j=1,i-1
1433           eps(i,j)=eps(j,i)
1434           epslip(i,j)=epslip(j,i)
1435         enddo
1436       enddo
1437       do i=1,ntyp
1438         do j=i,ntyp
1439           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1440           sigma(j,i)=sigma(i,j)
1441           rs0(i,j)=dwa16*sigma(i,j)
1442           rs0(j,i)=rs0(i,j)
1443         enddo
1444       enddo
1445       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1446      & 'Working parameters of the SC interactions:',
1447      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1448      & '  chi1   ','   chi2   '
1449       do i=1,ntyp
1450         do j=i,ntyp
1451           epsij=eps(i,j)
1452           epsijlip=epslip(i,j)
1453           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1454             rrij=sigma(i,j)
1455           else
1456             rrij=rr0(i)+rr0(j)
1457           endif
1458           r0(i,j)=rrij
1459           r0(j,i)=rrij
1460           rrij=rrij**expon
1461           epsij=eps(i,j)
1462           sigeps=dsign(1.0D0,epsij)
1463           epsij=dabs(epsij)
1464           aa_aq(i,j)=epsij*rrij*rrij
1465           bb_aq(i,j)=-sigeps*epsij*rrij
1466           aa_aq(j,i)=aa_aq(i,j)
1467           bb_aq(j,i)=bb_aq(i,j)
1468           sigeps=dsign(1.0D0,epsijlip)
1469           epsijlip=dabs(epsijlip)
1470           aa_lip(i,j)=epsijlip*rrij*rrij
1471           bb_lip(i,j)=-sigeps*epsijlip*rrij
1472           aa_lip(j,i)=aa_lip(i,j)
1473           bb_lip(j,i)=bb_lip(i,j)
1474           if (ipot.gt.2) then
1475             sigt1sq=sigma0(i)**2
1476             sigt2sq=sigma0(j)**2
1477             sigii1=sigii(i)
1478             sigii2=sigii(j)
1479             ratsig1=sigt2sq/sigt1sq
1480             ratsig2=1.0D0/ratsig1
1481             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1482             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1483             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1484           else
1485             rsum_max=sigma(i,j)
1486           endif
1487 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1488             sigmaii(i,j)=rsum_max
1489             sigmaii(j,i)=rsum_max
1490 c         else
1491 c           sigmaii(i,j)=r0(i,j)
1492 c           sigmaii(j,i)=r0(i,j)
1493 c         endif
1494 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1495           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1496             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1497             augm(i,j)=epsij*r_augm**(2*expon)
1498 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1499             augm(j,i)=augm(i,j)
1500           else
1501             augm(i,j)=0.0D0
1502             augm(j,i)=0.0D0
1503           endif
1504           if (lprint) then
1505             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1506      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1507      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1508           endif
1509         enddo
1510       enddo
1511 C
1512 C Define the SC-p interaction constants
1513 C
1514 #ifdef OLDSCP
1515       do i=1,20
1516 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1517 C helix formation)
1518 c       aad(i,1)=0.3D0*4.0D0**12
1519 C Following line for constants currently implemented
1520 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1521         aad(i,1)=1.5D0*4.0D0**12
1522 c       aad(i,1)=0.17D0*5.6D0**12
1523         aad(i,2)=aad(i,1)
1524 C "Soft" SC-p repulsion
1525         bad(i,1)=0.0D0
1526 C Following line for constants currently implemented
1527 c       aad(i,1)=0.3D0*4.0D0**6
1528 C "Hard" SC-p repulsion
1529         bad(i,1)=3.0D0*4.0D0**6
1530 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1531         bad(i,2)=bad(i,1)
1532 c       aad(i,1)=0.0D0
1533 c       aad(i,2)=0.0D0
1534 c       bad(i,1)=1228.8D0
1535 c       bad(i,2)=1228.8D0
1536       enddo
1537 #else
1538 C
1539 C 8/9/01 Read the SC-p interaction constants from file
1540 C
1541       do i=1,ntyp
1542         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1543       enddo
1544       do i=1,ntyp
1545         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1546         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1547         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1548         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1549       enddo
1550
1551       if (lprint) then
1552         write (iout,*) "Parameters of SC-p interactions:"
1553         do i=1,20
1554           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1555      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1556         enddo
1557       endif
1558 #endif
1559 C
1560 C Define the constants of the disulfide bridge
1561 C
1562 C      ebr=-12.0D0
1563 c
1564 c Old arbitrary potential - commented out.
1565 c
1566 c      dbr= 4.20D0
1567 c      fbr= 3.30D0
1568 c
1569 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1570 c energy surface of diethyl disulfide.
1571 c A. Liwo and U. Kozlowska, 11/24/03
1572 c
1573 C      D0CM = 3.78d0
1574 C      AKCM = 15.1d0
1575 C      AKTH = 11.0d0
1576 C      AKCT = 12.0d0
1577 C      V1SS =-1.08d0
1578 C      V2SS = 7.61d0
1579 C      V3SS = 13.7d0
1580       write (iout,*) dyn_ss,'dyndyn'
1581       if (dyn_ss) then
1582         ss_depth=ebr/wsc-0.25*eps(1,1)
1583 C        write(iout,*) akcm,whpb,wsc,'KURWA'
1584         Ht=Ht/wsc-0.25*eps(1,1)
1585
1586         akcm=akcm*whpb/wsc
1587         akth=akth*whpb/wsc
1588         akct=akct*whpb/wsc
1589         v1ss=v1ss*whpb/wsc
1590         v2ss=v2ss*whpb/wsc
1591         v3ss=v3ss*whpb/wsc
1592       else
1593         ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1594       endif
1595
1596 C      if (lprint) then
1597       write (iout,'(/a)') "Disulfide bridge parameters:"
1598       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1599       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1600       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1601       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1602      & ' v3ss:',v3ss
1603 C      endif
1604       if (shield_mode.gt.0) then
1605       pi=3.141592d0
1606 C VSolvSphere the volume of solving sphere
1607 C      print *,pi,"pi"
1608 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1609 C there will be no distinction between proline peptide group and normal peptide
1610 C group in case of shielding parameters
1611       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1612       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1613       write (iout,*) "VSolvSphere",VSolvSphere," VSolvSphere_div",
1614      &  VSolvSphere_div
1615 C long axis of side chain
1616       do i=1,ntyp
1617       long_r_sidechain(i)=vbldsc0(1,i)
1618       short_r_sidechain(i)=sigma0(i)
1619       enddo
1620       buff_shield=1.0d0
1621       endif
1622       return
1623   111 write (iout,*) "Error reading bending energy parameters."
1624       goto 999
1625   112 write (iout,*) "Error reading rotamer energy parameters."
1626       goto 999
1627   113 write (iout,*) "Error reading torsional energy parameters."
1628       goto 999
1629   114 write (iout,*) "Error reading double torsional energy parameters."
1630       goto 999
1631   115 write (iout,*)
1632      &  "Error reading cumulant (multibody energy) parameters."
1633       goto 999
1634   116 write (iout,*) "Error reading electrostatic energy parameters."
1635       goto 999
1636  1161 write (iout,*) "Error reading lipid parameters."
1637       goto 999
1638   117 write (iout,*) "Error reading side chain interaction parameters."
1639       goto 999
1640   118 write (iout,*) "Error reading SCp interaction parameters."
1641       goto 999
1642   119 write (iout,*) "Error reading SCCOR parameters"
1643       goto 999
1644   121 write (iout,*) "Error reading bond parameters"
1645   999 continue
1646 #ifdef MPI
1647       call MPI_Finalize(Ierror)
1648 #endif
1649       stop
1650       return
1651       end