source/maxlik/src_MD_T_maxlik-NEWCORR/parmread.F

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