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