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