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