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