source/maxlik/src-Fmatch_safe/parmread.F.safe

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