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       include 'COMMON.CONTROL'
  25       character*1 t1,t2,t3
  26       character*1 onelett(4) /"G","A","P","D"/
  27       character*1 toronelet(-2:2) /"p","a","G","A","P"/
  28       logical lprint
  29       dimension blower(3,3,maxlob)
  30       character*800 controlcard
  31       character*256 bondname_t,thetname_t,rotname_t,torname_t,
  32      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
  33      & sccorname_t
  34       integer ilen
  35       external ilen
  36       character*16 key
  37       integer iparm
  38       double precision ip,mp
  39       character*6 res1
  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
 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,vbldpdum,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,vbldpdum,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       read (ithep,*) nthetyp,ntheterm,ntheterm2,
 361      &  ntheterm3,nsingle,ndouble
 362       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 363       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 364       do i=-ntyp1,-1
 365         ithetyp(i)=-ithetyp(-i)
 366       enddo
 367 c      write (iout,*) "tu dochodze"
 368       do iblock=1,2
 369       do i=-maxthetyp,maxthetyp
 370         do j=-maxthetyp,maxthetyp
 371           do k=-maxthetyp,maxthetyp
 372             aa0thet(i,j,k,iblock)=0.0d0
 373             do l=1,ntheterm
 374               aathet(l,i,j,k,iblock)=0.0d0
 375             enddo
 376             do l=1,ntheterm2
 377               do m=1,nsingle
 378                 bbthet(m,l,i,j,k,iblock)=0.0d0
 379                 ccthet(m,l,i,j,k,iblock)=0.0d0
 380                 ddthet(m,l,i,j,k,iblock)=0.0d0
 381                 eethet(m,l,i,j,k,iblock)=0.0d0
 382               enddo
 383             enddo
 384             do l=1,ntheterm3
 385               do m=1,ndouble
 386                 do mm=1,ndouble
 387                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 388                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 389                 enddo
 390               enddo
 391             enddo
 392           enddo
 393         enddo
 394       enddo
 395       enddo
 396 C      write (iout,*) "KURWA1"
 397       do iblock=1,2
 398       do i=0,nthetyp
 399         do j=-nthetyp,nthetyp
 400           do k=-nthetyp,nthetyp
 401             read (ithep,'(6a)') res1
 402 c            write(iout,*) res1,i,j,k
 403             read (ithep,*) aa0thet(i,j,k,iblock)
 404             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 405             read (ithep,*)
 406      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 407      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 408      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 409      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 410      &        ,ll=1,ntheterm2)
 411             read (ithep,*)
 412      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 413      &      ffthet(lll,llll,ll,i,j,k,iblock),
 414      &         ggthet(llll,lll,ll,i,j,k,iblock)
 415      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 416      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 417           enddo
 418         enddo
 419       enddo
 420 C
 421 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 422 C coefficients of theta-and-gamma-dependent terms are zero.
 423 C
 424       do i=1,nthetyp
 425         do j=1,nthetyp
 426           do l=1,ntheterm
 427             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 428             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 429           enddo
 430           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 431           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 432         enddo
 433         do l=1,ntheterm
 434           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 435         enddo
 436         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 437       enddo
 438       enddo
 439 C       write(iout,*) "KURWA1.5"
 440 C Substitution for D aminoacids from symmetry.
 441       do iblock=1,2
 442       do i=-nthetyp,0
 443         do j=-nthetyp,nthetyp
 444           do k=-nthetyp,nthetyp
 445            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 446            do l=1,ntheterm
 447            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 448            enddo
 449            do ll=1,ntheterm2
 450             do lll=1,nsingle
 451             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 452             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 453             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 454             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 455             enddo
 456           enddo
 457           do ll=1,ntheterm3
 458            do lll=2,ndouble
 459             do llll=1,lll-1
 460             ffthet(llll,lll,ll,i,j,k,iblock)=
 461      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 462             ffthet(lll,llll,ll,i,j,k,iblock)=
 463      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 464             ggthet(llll,lll,ll,i,j,k,iblock)=
 465      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 466             ggthet(lll,llll,ll,i,j,k,iblock)=
 467      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 468             enddo !ll
 469            enddo  !lll
 470           enddo   !llll
 471          enddo    !k
 472         enddo     !j
 473        enddo      !i
 474       enddo       !iblock
 475
 476 C
 477 C Control printout of the coefficients of virtual-bond-angle potentials
 478 C
 479       if (lprint) then
 480         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 481         do i=1,nthetyp+1
 482           do j=1,nthetyp+1
 483             do k=1,nthetyp+1
 484               write (iout,'(//4a)')
 485      &         'Type ',onelett(i),onelett(j),onelett(k)
 486               write (iout,'(//a,10x,a)') " l","a[l]"
 487               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 488               write (iout,'(i2,1pe15.5)')
 489      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 490             do l=1,ntheterm2
 491               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 492      &          "b",l,"c",l,"d",l,"e",l
 493               do m=1,nsingle
 494                 write (iout,'(i2,4(1pe15.5))') m,
 495      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 496      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 497               enddo
 498             enddo
 499             do l=1,ntheterm3
 500               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 501      &          "f+",l,"f-",l,"g+",l,"g-",l
 502               do m=2,ndouble
 503                 do n=1,m-1
 504                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 505      &              ffthet(n,m,l,i,j,k,iblock),
 506      &              ffthet(m,n,l,i,j,k,iblock),
 507      &              ggthet(n,m,l,i,j,k,iblock),
 508      &              ggthet(m,n,l,i,j,k,iblock)
 509                 enddo
 510               enddo
 511             enddo
 512           enddo
 513         enddo
 514       enddo
 515       call flush(iout)
 516       endif
 517 #endif
 518
 519 #ifdef CRYST_SC
 520 C
 521 C Read the parameters of the probability distribution/energy expression
 522 C of the side chains.
 523 C
 524       do i=1,ntyp
 525         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 526         if (i.eq.10) then
 527           dsc_inv(i)=0.0D0
 528         else
 529           dsc_inv(i)=1.0D0/dsc(i)
 530         endif
 531         if (i.ne.10) then
 532         do j=1,nlob(i)
 533           do k=1,3
 534             do l=1,3
 535               blower(l,k,j)=0.0D0
 536             enddo
 537           enddo
 538         enddo
 539         bsc(1,i)=0.0D0
 540         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 541         censc(1,1,-i)=censc(1,1,i)
 542         censc(2,1,-i)=censc(2,1,i)
 543         censc(3,1,-i)=-censc(3,1,i)
 544         do j=2,nlob(i)
 545           read (irotam,*) bsc(j,i)
 546           read (irotam,*) (censc(k,j,i),k=1,3),
 547      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 548         censc(1,j,-i)=censc(1,j,i)
 549         censc(2,j,-i)=censc(2,j,i)
 550         censc(3,j,-i)=-censc(3,j,i)
 551 C BSC is amplitude of Gaussian
 552         enddo
 553         do j=1,nlob(i)
 554           do k=1,3
 555             do l=1,k
 556               akl=0.0D0
 557               do m=1,3
 558                 akl=akl+blower(k,m,j)*blower(l,m,j)
 559               enddo
 560               gaussc(k,l,j,i)=akl
 561               gaussc(l,k,j,i)=akl
 562              if (((k.eq.3).and.(l.ne.3))
 563      &        .or.((l.eq.3).and.(k.ne.3))) then
 564                 gaussc(k,l,j,-i)=-akl
 565                 gaussc(l,k,j,-i)=-akl
 566               else
 567                 gaussc(k,l,j,-i)=akl
 568                 gaussc(l,k,j,-i)=akl
 569               endif
 570             enddo
 571           enddo
 572         enddo
 573         endif
 574       enddo
 575       close (irotam)
 576       if (lprint) then
 577         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 578         do i=1,ntyp
 579           nlobi=nlob(i)
 580           if (nlobi.gt.0) then
 581           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 582      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 583 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 584 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 585 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 586 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 587            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 588      &                             'log h',(bsc(j,i),j=1,nlobi)
 589            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 590      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 591 c          write (iout,'(a)')
 592 c         do j=1,nlobi
 593 c           ind=0
 594 c           do k=1,3
 595 c             do l=1,k
 596 c              ind=ind+1
 597 c              blower(k,l,j)=gaussc(ind,j,i)
 598 c             enddo
 599 c           enddo
 600 c         enddo
 601           do k=1,3
 602             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 603      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 604           enddo
 605           endif
 606         enddo
 607       endif
 608 #else
 609 C
 610 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 611 C added by Urszula Kozlowska 07/11/2007
 612 C
 613       do i=1,ntyp
 614         read (irotam,*)
 615        if (i.eq.10) then
 616          read (irotam,*)
 617        else
 618          do j=1,65
 619            read(irotam,*) sc_parmin(j,i)
 620          enddo
 621        endif
 622       enddo
 623 #endif
 624       close(irotam)
 625 #ifdef CRYST_TOR
 626 C
 627 C Read torsional parameters in old format
 628 C
 629       read (itorp,*) ntortyp,nterm_old
 630       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 631       read (itorp,*) (itortyp(i),i=1,ntyp)
 632       do i=1,ntortyp
 633         do j=1,ntortyp
 634           read (itorp,'(a)')
 635           do k=1,nterm_old
 636             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 637           enddo
 638         enddo
 639       enddo
 640       close (itorp)
 641       if (lprint) then
 642         write (iout,'(/a/)') 'Torsional constants:'
 643         do i=1,ntortyp
 644           do j=1,ntortyp
 645             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 646             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 647           enddo
 648         enddo
 649       endif
 650
 651
 652 #else
 653 C
 654 C Read torsional parameters
 655 C
 656       read (itorp,*) ntortyp
 657       read (itorp,*) (itortyp(i),i=1,ntyp)
 658       write (iout,*) 'ntortyp',ntortyp
 659       do iblock=1,2
 660       do i=-ntyp,-1
 661        itortyp(i)=-itortyp(-i)
 662       enddo
 663 c      write (iout,*) 'ntortyp',ntortyp
 664       do i=0,ntortyp-1
 665         do j=-ntortyp+1,ntortyp-1
 666           read (itorp,*) nterm(i,j,iblock),
 667      &          nlor(i,j,iblock)
 668           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 669           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 670           v0ij=0.0d0
 671           si=-1.0d0
 672           do k=1,nterm(i,j,iblock)
 673             read (itorp,*) kk,v1(k,i,j,iblock),
 674      &      v2(k,i,j,iblock)
 675             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 676             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 677             v0ij=v0ij+si*v1(k,i,j,iblock)
 678             si=-si
 679          enddo
 680           do k=1,nlor(i,j,iblock)
 681             read (itorp,*) kk,vlor1(k,i,j),
 682      &        vlor2(k,i,j),vlor3(k,i,j)
 683             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 684           enddo
 685           v0(i,j,iblock)=v0ij
 686           v0(-i,-j,iblock)=v0ij
 687         enddo
 688       enddo
 689       enddo
 690       close (itorp)
 691       if (lprint) then
 692         write (iout,'(/a/)') 'Torsional constants:'
 693         do i=1,ntortyp
 694           do j=1,ntortyp
 695             write (iout,*) 'ityp',i,' jtyp',j
 696             write (iout,*) 'Fourier constants'
 697             do k=1,nterm(i,j,iblock)
 698               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 699      &        v2(k,i,j,iblock)
 700             enddo
 701             write (iout,*) 'Lorenz constants'
 702             do k=1,nlor(i,j,iblock)
 703               write (iout,'(3(1pe15.5))')
 704      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 705             enddo
 706           enddo
 707         enddo
 708       endif
 709 C
 710 C 6/23/01 Read parameters for double torsionals
 711 C
 712       do iblock=1,2
 713       do i=0,ntortyp-1
 714         do j=-ntortyp+1,ntortyp-1
 715           do k=-ntortyp+1,ntortyp-1
 716             read (itordp,'(3a1)') t1,t2,t3
 717 c              write (iout,*) "OK onelett",
 718 c     &         i,j,k,t1,t2,t3
 719
 720             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 721      &        .or. t3.ne.toronelet(k)) then
 722               write (iout,*) "Error in double torsional parameter file",
 723      &         i,j,k,t1,t2,t3
 724 #ifdef MPI
 725               call MPI_Finalize(Ierror)
 726 #endif
 727                stop "Error in double torsional parameter file"
 728             endif
 729           read (itordp,*) ntermd_1(i,j,k,iblock),
 730      &         ntermd_2(i,j,k,iblock)
 731             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 732             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 733             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
 734      &         ntermd_1(i,j,k,iblock))
 735             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
 736      &         ntermd_1(i,j,k,iblock))
 737             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
 738      &         ntermd_1(i,j,k,iblock))
 739             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
 740      &         ntermd_1(i,j,k,iblock))
 741 C Martix of D parameters for one dimesional foureir series
 742             do l=1,ntermd_1(i,j,k,iblock)
 743              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 744              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 745              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 746              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 747 c            write(iout,*) "whcodze" ,
 748 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 749             enddo
 750             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
 751      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 752      &         v2s(m,l,i,j,k,iblock),
 753      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 754 C Martix of D parameters for two dimesional fourier series
 755             do l=1,ntermd_2(i,j,k,iblock)
 756              do m=1,l-1
 757              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 758              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 759              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 760              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 761              enddo!m
 762             enddo!l
 763           enddo!k
 764         enddo!j
 765       enddo!i
 766       enddo!iblock
 767       if (lprint) then
 768       write (iout,*)
 769       write (iout,*) 'Constants for double torsionals'
 770       do iblock=1,2
 771       do i=0,ntortyp-1
 772         do j=-ntortyp+1,ntortyp-1
 773           do k=-ntortyp+1,ntortyp-1
 774             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 775      &        ' nsingle',ntermd_1(i,j,k,iblock),
 776      &        ' ndouble',ntermd_2(i,j,k,iblock)
 777             write (iout,*)
 778             write (iout,*) 'Single angles:'
 779             do l=1,ntermd_1(i,j,k,iblock)
 780               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 781      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 782      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 783      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 784             enddo
 785             write (iout,*)
 786             write (iout,*) 'Pairs of angles:'
 787             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 788             do l=1,ntermd_2(i,j,k,iblock)
 789               write (iout,'(i5,20f10.5)')
 790      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 791             enddo
 792             write (iout,*)
 793            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 794             do l=1,ntermd_2(i,j,k,iblock)
 795               write (iout,'(i5,20f10.5)')
 796      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 797      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 798             enddo
 799             write (iout,*)
 800           enddo
 801         enddo
 802       enddo
 803       enddo
 804       endif
 805 #endif
 806 C Read of Side-chain backbone correlation parameters
 807 C Modified 11 May 2012 by Adasko
 808 CCC
 809 C
 810       read (isccor,*) nsccortyp
 811       read (isccor,*) (isccortyp(i),i=1,ntyp)
 812       do i=-ntyp,-1
 813         isccortyp(i)=-isccortyp(-i)
 814       enddo
 815       iscprol=isccortyp(20)
 816 c      write (iout,*) 'ntortyp',ntortyp
 817       maxinter=3
 818 cc maxinter is maximum interaction sites
 819       do l=1,maxinter
 820       do i=1,nsccortyp
 821         do j=1,nsccortyp
 822           read (isccor,*)
 823      &nterm_sccor(i,j),nlor_sccor(i,j)
 824 c          write (iout,*) nterm_sccor(i,j)
 825           v0ijsccor=0.0d0
 826           v0ijsccor1=0.0d0
 827           v0ijsccor2=0.0d0
 828           v0ijsccor3=0.0d0
 829           si=-1.0d0
 830           nterm_sccor(-i,j)=nterm_sccor(i,j)
 831           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 832           nterm_sccor(i,-j)=nterm_sccor(i,j)
 833 c          write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
 834 c     &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
 835           do k=1,nterm_sccor(i,j)
 836             read (isccor,*) kk,v1sccor(k,l,i,j)
 837      &    ,v2sccor(k,l,i,j)
 838             if (j.eq.iscprol) then
 839              if (i.eq.isccortyp(10)) then
 840              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 841              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 842              else
 843              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 844      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 845              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 846      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 847              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 848              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 849              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 850              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 851              endif
 852             else
 853              if (i.eq.isccortyp(10)) then
 854              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 855              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 856              else
 857                if (j.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              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 862              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 863              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 864              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 865              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 866              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 867                 endif
 868                endif
 869             endif
 870             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 871             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 872             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 873             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 874             si=-si
 875            enddo
 876           do k=1,nlor_sccor(i,j)
 877             read (isccor,*) kk,vlor1sccor(k,i,j),
 878      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 879             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 880      &(1+vlor3sccor(k,i,j)**2)
 881           enddo
 882           v0sccor(l,i,j)=v0ijsccor
 883           v0sccor(l,-i,j)=v0ijsccor1
 884           v0sccor(l,i,-j)=v0ijsccor2
 885           v0sccor(l,-i,-j)=v0ijsccor3
 886           enddo
 887         enddo
 888       enddo
 889       close (isccor)
 890       if (lprint) then
 891         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 892         do i=1,nsccortyp
 893           do j=1,nsccortyp
 894             write (iout,*) 'ityp',i,' jtyp',j
 895             write (iout,*) 'Fourier constants'
 896             do k=1,nterm_sccor(i,j)
 897               write (iout,'(2(1pe15.5))')
 898      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 899             enddo
 900             write (iout,*) 'Lorenz constants'
 901             do k=1,nlor_sccor(i,j)
 902               write (iout,'(3(1pe15.5))')
 903      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 904             enddo
 905           enddo
 906         enddo
 907       endif
 908 C
 909 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 910 C         interaction energy of the Gly, Ala, and Pro prototypes.
 911 C
 912       read (ifourier,*) nloctyp
 913       do i=0,nloctyp-1
 914         read (ifourier,*)
 915         read (ifourier,*) (b(ii,i),ii=1,13)
 916         if (lprint) then
 917         write (iout,*) 'Type',i
 918         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 919         endif
 920         B1(1,i)  = b(3,i)
 921         B1(2,i)  = b(5,i)
 922         B1(1,-i) = b(3,i)
 923         B1(2,-i) = -b(5,i)
 924 c        b1(1,i)=0.0d0
 925 c        b1(2,i)=0.0d0
 926         B1tilde(1,i) = b(3,i)
 927         B1tilde(2,i) =-b(5,i)
 928         B1tilde(1,-i) =-b(3,i)
 929         B1tilde(2,-i) =b(5,i)
 930 c        b1tilde(1,i)=0.0d0
 931 c        b1tilde(2,i)=0.0d0
 932         B2(1,i)  = b(2,i)
 933         B2(2,i)  = b(4,i)
 934         B2(1,-i)  =b(2,i)
 935         B2(2,-i)  =-b(4,i)
 936
 937 c        b2(1,i)=0.0d0
 938 c        b2(2,i)=0.0d0
 939         CC(1,1,i)= b(7,i)
 940         CC(2,2,i)=-b(7,i)
 941         CC(2,1,i)= b(9,i)
 942         CC(1,2,i)= b(9,i)
 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 c        CC(1,1,i)=0.0d0
 948 c        CC(2,2,i)=0.0d0
 949 c        CC(2,1,i)=0.0d0
 950 c        CC(1,2,i)=0.0d0
 951         Ctilde(1,1,i)=b(7,i)
 952         Ctilde(1,2,i)=b(9,i)
 953         Ctilde(2,1,i)=-b(9,i)
 954         Ctilde(2,2,i)=b(7,i)
 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
 960 c        Ctilde(1,1,i)=0.0d0
 961 c        Ctilde(1,2,i)=0.0d0
 962 c        Ctilde(2,1,i)=0.0d0
 963 c        Ctilde(2,2,i)=0.0d0
 964         DD(1,1,i)= b(6,i)
 965         DD(2,2,i)=-b(6,i)
 966         DD(2,1,i)= b(8,i)
 967         DD(1,2,i)= b(8,i)
 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 c        DD(1,1,i)=0.0d0
 973 c        DD(2,2,i)=0.0d0
 974 c        DD(2,1,i)=0.0d0
 975 c        DD(1,2,i)=0.0d0
 976         Dtilde(1,1,i)=b(6,i)
 977         Dtilde(1,2,i)=b(8,i)
 978         Dtilde(2,1,i)=-b(8,i)
 979         Dtilde(2,2,i)=b(6,i)
 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
 985 c        Dtilde(1,1,i)=0.0d0
 986 c        Dtilde(1,2,i)=0.0d0
 987 c        Dtilde(2,1,i)=0.0d0
 988 c        Dtilde(2,2,i)=0.0d0
 989         EE(1,1,i)= b(10,i)+b(11,i)
 990         EE(2,2,i)=-b(10,i)+b(11,i)
 991         EE(2,1,i)= b(12,i)-b(13,i)
 992         EE(1,2,i)= b(12,i)+b(13,i)
 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
 998 c        ee(1,1,i)=1.0d0
 999 c        ee(2,2,i)=1.0d0
1000 c        ee(2,1,i)=0.0d0
1001 c        ee(1,2,i)=0.0d0
1002 c        ee(2,1,i)=ee(1,2,i)
1003
1004       enddo
1005       if (lprint) then
1006       do i=1,nloctyp
1007         write (iout,*) 'Type',i
1008         write (iout,*) 'B1'
1009 c        write (iout,'(f10.5)') B1(:,i)
1010         write(iout,*) B1(1,i),B1(2,i)
1011         write (iout,*) 'B2'
1012 c        write (iout,'(f10.5)') B2(:,i)
1013         write(iout,*) B2(1,i),B2(2,i)
1014         write (iout,*) 'CC'
1015         do j=1,2
1016           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1017         enddo
1018         write(iout,*) 'DD'
1019         do j=1,2
1020           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1021         enddo
1022         write(iout,*) 'EE'
1023         do j=1,2
1024           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1025         enddo
1026       enddo
1027       endif
1028 C
1029 C Read electrostatic-interaction parameters
1030 C
1031       if (lprint) then
1032         write (iout,'(/a)') 'Electrostatic interaction constants:'
1033         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1034      &            'IT','JT','APP','BPP','AEL6','AEL3'
1035       endif
1036       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1037       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1038       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1039       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1040       close (ielep)
1041       do i=1,2
1042         do j=1,2
1043         rri=rpp(i,j)**6
1044         app (i,j)=epp(i,j)*rri*rri
1045         bpp (i,j)=-2.0D0*epp(i,j)*rri
1046         ael6(i,j)=elpp6(i,j)*4.2D0**6
1047         ael3(i,j)=elpp3(i,j)*4.2D0**3
1048         lprint=.true.
1049         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1050      &                    ael6(i,j),ael3(i,j)
1051          lprint=.false.
1052         enddo
1053       enddo
1054 C
1055 C Read side-chain interaction parameters.
1056 C
1057       read (isidep,*) ipot,expon
1058       if (ipot.lt.1 .or. ipot.gt.5) then
1059         write (iout,'(2a)') 'Error while reading SC interaction',
1060      &               'potential file - unknown potential type.'
1061         stop
1062       endif
1063       expon2=expon/2
1064       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1065      & ', exponents are ',expon,2*expon
1066       goto (10,20,30,30,40) ipot
1067 C----------------------- LJ potential ---------------------------------
1068    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1069       if (lprint) then
1070         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1071         write (iout,'(a/)') 'The epsilon array:'
1072         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1073         write (iout,'(/a)') 'One-body parameters:'
1074         write (iout,'(a,4x,a)') 'residue','sigma'
1075         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1076       endif
1077       goto 50
1078 C----------------------- LJK potential --------------------------------
1079    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1080      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1081       if (lprint) then
1082         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1083         write (iout,'(a/)') 'The epsilon array:'
1084         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1085         write (iout,'(/a)') 'One-body parameters:'
1086         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1087         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1088      &        i=1,ntyp)
1089       endif
1090       goto 50
1091 C---------------------- GB or BP potential -----------------------------
1092    30 do i=1,ntyp
1093        read (isidep,*)(eps(i,j),j=i,ntyp)
1094       enddo
1095       read (isidep,*)(sigma0(i),i=1,ntyp)
1096       read (isidep,*)(sigii(i),i=1,ntyp)
1097       read (isidep,*)(chip(i),i=1,ntyp)
1098       read (isidep,*)(alp(i),i=1,ntyp)
1099       do i=1,ntyp
1100        read (isidep,*)(epslip(i,j),j=i,ntyp)
1101 C       write(iout,*) "WARNING!!",i,ntyp
1102        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1103 C       do j=1,ntyp
1104 C       epslip(i,j)=epslip(i,j)+0.05d0
1105 C       enddo
1106       enddo
1107 C For the GB potential convert sigma'**2 into chi'
1108       if (ipot.eq.4) then
1109         do i=1,ntyp
1110           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1111         enddo
1112       endif
1113       if (lprint) then
1114         write (iout,'(/a/)') 'Parameters of the BP potential:'
1115         write (iout,'(a/)') 'The epsilon array:'
1116         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1117         write (iout,'(/a)') 'One-body parameters:'
1118         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1119      &       '    chip  ','    alph  '
1120         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1121      &                     chip(i),alp(i),i=1,ntyp)
1122       endif
1123       goto 50
1124 C--------------------- GBV potential -----------------------------------
1125    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1126      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1127      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1128       if (lprint) then
1129         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1130         write (iout,'(a/)') 'The epsilon array:'
1131         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1132         write (iout,'(/a)') 'One-body parameters:'
1133         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1134      &      's||/s_|_^2','    chip  ','    alph  '
1135         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1136      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1137       endif
1138    50 continue
1139       close (isidep)
1140 C-----------------------------------------------------------------------
1141 C Calculate the "working" parameters of SC interactions.
1142       do i=2,ntyp
1143         do j=1,i-1
1144           eps(i,j)=eps(j,i)
1145           epslip(i,j)=epslip(j,i)
1146         enddo
1147       enddo
1148       do i=1,ntyp
1149         do j=i,ntyp
1150           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1151           sigma(j,i)=sigma(i,j)
1152           rs0(i,j)=dwa16*sigma(i,j)
1153           rs0(j,i)=rs0(i,j)
1154         enddo
1155       enddo
1156       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1157      & 'Working parameters of the SC interactions:',
1158      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1159      & '  chi1   ','   chi2   '
1160       do i=1,ntyp
1161         do j=i,ntyp
1162           epsij=eps(i,j)
1163           epsijlip=epslip(i,j)
1164           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1165             rrij=sigma(i,j)
1166           else
1167             rrij=rr0(i)+rr0(j)
1168           endif
1169           r0(i,j)=rrij
1170           r0(j,i)=rrij
1171           rrij=rrij**expon
1172           epsij=eps(i,j)
1173           sigeps=dsign(1.0D0,epsij)
1174           epsij=dabs(epsij)
1175           aa_aq(i,j)=epsij*rrij*rrij
1176           bb_aq(i,j)=-sigeps*epsij*rrij
1177           aa_aq(j,i)=aa_aq(i,j)
1178           bb_aq(j,i)=bb_aq(i,j)
1179           sigeps=dsign(1.0D0,epsijlip)
1180           epsijlip=dabs(epsijlip)
1181           aa_lip(i,j)=epsijlip*rrij*rrij
1182           bb_lip(i,j)=-sigeps*epsijlip*rrij
1183           aa_lip(j,i)=aa_lip(i,j)
1184           bb_lip(j,i)=bb_lip(i,j)
1185           if (ipot.gt.2) then
1186             sigt1sq=sigma0(i)**2
1187             sigt2sq=sigma0(j)**2
1188             sigii1=sigii(i)
1189             sigii2=sigii(j)
1190             ratsig1=sigt2sq/sigt1sq
1191             ratsig2=1.0D0/ratsig1
1192             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1193             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1194             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1195           else
1196             rsum_max=sigma(i,j)
1197           endif
1198 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1199             sigmaii(i,j)=rsum_max
1200             sigmaii(j,i)=rsum_max
1201 c         else
1202 c           sigmaii(i,j)=r0(i,j)
1203 c           sigmaii(j,i)=r0(i,j)
1204 c         endif
1205 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1206           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1207             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1208             augm(i,j)=epsij*r_augm**(2*expon)
1209 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1210             augm(j,i)=augm(i,j)
1211           else
1212             augm(i,j)=0.0D0
1213             augm(j,i)=0.0D0
1214           endif
1215           if (lprint) then
1216             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1217      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1218      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1219           endif
1220         enddo
1221       enddo
1222 C
1223 C Define the SC-p interaction constants
1224 C
1225 #ifdef OLDSCP
1226       do i=1,20
1227 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1228 C helix formation)
1229 c       aad(i,1)=0.3D0*4.0D0**12
1230 C Following line for constants currently implemented
1231 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1232         aad(i,1)=1.5D0*4.0D0**12
1233 c       aad(i,1)=0.17D0*5.6D0**12
1234         aad(i,2)=aad(i,1)
1235 C "Soft" SC-p repulsion
1236         bad(i,1)=0.0D0
1237 C Following line for constants currently implemented
1238 c       aad(i,1)=0.3D0*4.0D0**6
1239 C "Hard" SC-p repulsion
1240         bad(i,1)=3.0D0*4.0D0**6
1241 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1242         bad(i,2)=bad(i,1)
1243 c       aad(i,1)=0.0D0
1244 c       aad(i,2)=0.0D0
1245 c       bad(i,1)=1228.8D0
1246 c       bad(i,2)=1228.8D0
1247       enddo
1248 #else
1249 C
1250 C 8/9/01 Read the SC-p interaction constants from file
1251 C
1252       do i=1,ntyp
1253         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1254       enddo
1255       do i=1,ntyp
1256         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1257         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1258         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1259         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1260       enddo
1261
1262       if (lprint) then
1263         write (iout,*) "Parameters of SC-p interactions:"
1264         do i=1,20
1265           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1266      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1267         enddo
1268       endif
1269 #endif
1270 C
1271 C Define the constants of the disulfide bridge
1272 C
1273 C      ebr=-12.0D0
1274 c
1275 c Old arbitrary potential - commented out.
1276 c
1277 c      dbr= 4.20D0
1278 c      fbr= 3.30D0
1279 c
1280 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1281 c energy surface of diethyl disulfide.
1282 c A. Liwo and U. Kozlowska, 11/24/03
1283 c
1284 C      D0CM = 3.78d0
1285 C      AKCM = 15.1d0
1286 C      AKTH = 11.0d0
1287 C      AKCT = 12.0d0
1288 C      V1SS =-1.08d0
1289 C      V2SS = 7.61d0
1290 C      V3SS = 13.7d0
1291       write (iout,*) dyn_ss,'dyndyn'
1292       if (dyn_ss) then
1293         ss_depth=ebr/wsc-0.25*eps(1,1)
1294 C        write(iout,*) akcm,whpb,wsc,'KURWA'
1295         Ht=Ht/wsc-0.25*eps(1,1)
1296
1297         akcm=akcm*whpb/wsc
1298         akth=akth*whpb/wsc
1299         akct=akct*whpb/wsc
1300         v1ss=v1ss*whpb/wsc
1301         v2ss=v2ss*whpb/wsc
1302         v3ss=v3ss*whpb/wsc
1303       else
1304         ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1305       endif
1306
1307 C      if (lprint) then
1308       write (iout,'(/a)') "Disulfide bridge parameters:"
1309       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1310       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1311       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1312       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1313      & ' v3ss:',v3ss
1314 C      endif
1315       return
1316       end